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Executive  Summary 


This  work  performed  a  review  of  wide-band  DOA  algorithms  in  the  literature 
which  was  accumulated  for  more  than  30  years.  There  are  more  than  fifty  publications 
for  wide-band  detection  of  Direction  of  Arrival  (DOA)  algorithms  which  are  available  in 
the  literature.  We  have  reviewed  the  most  relevant  one  and  have  not  reviewed  others 
which  will  not  be  applicable  or  suitable  for  hardware  implementation.  These  algorithms 
were  generally  presented  in  a  very  complex  or  condense  form,  which  are  not  easily 
understandable  for  people  who  are  outside  that  narrow  field.  One  of  the  reasons  for  their 
complex  representation  is  due  to  their  publication  in  IEEE  transactions  and  conferences. 
These  transactions  generally  prefer  highly  mathematical  papers  and  sometimes  authors 
insert  mathematics  so  chances  of  their  papers  are  increased.  Another  reason  for  condense 
reporting  is  that  these  papers  face  a  page  limit.  Therefore  algorithms  need  to  be 
accommodated  within  those  guidelines  and  also  comply  with  the  reviewers  comments. 
One  unfortunate  thing  happens  in  this  process  that  essential  information  does  not  get  into 
the  papers  and  there  is  always  missing  information.  This  missing  information  is  acute  in 
our  case  as  we  are  looking  from  hardware  implementation  point  of  view  and  we  are 
ignoring  details  of  statistical  results  and  errors  which  are  irrelevant  in  our  case.  We  are 
willing  to  sacrifice  small  amount  of  error  in  order  to  accomplish  the  goal  of  implementing 
them  in  hardware  for  real  time  applications. 
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We  have  discovered  a  class  of  computational  requirements  that  would  be  required 
in  all  these  algorithms  as  was  summarized.  We  have  also  given  reviews  and  challenges. 
We  were  able  to  cut  through  all  the  mathematics  and  convert  algorithms  into  simple 
arithmetic  operations.  This  step  is  very  useful  in  visualizing  an  architecture.  We  have 
filled  a  gap  between  the  design  of  computer  hardware  especially  special  purpose  parallel 
architectures  and  available  algorithms  for  various  wide-band  DOA  algorithms. 

This  work  was  the  first  step  in  sorting  out  which  algorithm  is  appropriate  for 
further  study  and  for  its  hardware  implementation  for  real  time  applications.  We  have 
developed  hardware  as  described  in  Chapter  5  and  Volume  2  of  this  report.  Work  is  in 
progress  for  implementation  of  identified  computational  steps  in  FPGAs. 

This  work  can  be  extended  to  develop  re-configurable  test-bed  environment  for 
investigative  studies  for  various  algorithm.  The  re-configurable  test-bed  would  be  useful 
to  study  timing,  memory,  hardware  requirement  and  accuracy  of  results  for  various 
algorithms.  This  test-bed  would  be  useful  in  evaluating  different  number  of  sensors  and 
different  kind  of  sensors.  This  test-bed  would  also  be  a  technology  scalable  system  and 
would  become  useful  in  deployment  hardware. 
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Chapter  1 


Introduction 

Array  processing  has  been  an  important  part  of  signal  processing  in  the  past  few 
decades  [1-40].  The  array  consists  of  sensors  located  at  different  points  in  space  with 
respect  to  a  reference  point.  Direction  of  Arrival  (DOA)  denotes  the  direction  from  which 
the  wave  fields  arrive  at  the  sensor  array.  The  goal  in  DOA  detection  and  estimation  is  to 
accurately  determine  the  number  of  sources  producing  waveforms  and  the  locations  of 
those  sources.  The  passive  detection  of  objects  has  become  important  in  the  military 
applications  as  it  evades  detection  by  others.  The  estimation  of  Direction-Of-Arrival 
(DOA)  from  energy  wavefield  has  many  applications  both  inside  and  outside  of  the 
military  use.  Some  civilian  applications  are  in  the  areas  of  communications,  air  traffic 
control,  seismology,  sonar,  and  bioengineering.  Generally  passive  detection  approaches 
are  computed  intensive.  Their  hardware  implementation  could  be  cumbersome  and  hard 
to  meet  the  real  time  computational  requirements. 

With  the  growth  in  technology  and  increase  in  processing  power,  it  is  now 
possible  to  develop  real  time  hardware  for  many  applications.  There  are  three  types  of 
techniques  generally  available  in  the  literature  for  detection  of  DOA  [1-6].  They  are 
Power  Spectral  (PSD),  maximum  likelihood,  maximum  entropy  and  signal  subspace 
methods.  PSD  techniques  are  ineffective  and  maximum  likelihood  techniques  are 
considered  accurate,  they  are  very  computationally  expensive  and  also  yield  non-linear 
equations.  Maximum  entropy  methods  introduce  bias  and  are  very  sensitive  to  various 
parameters.  Signal  subspace  techniques  are  accurate,  give  high  resolution  and  reasonably 
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compute  intensive.  Subspace-based  detection  techniques  that  yield  high  resolution  for 
estimating  DOA  through  passive  sensors  have  been  popular  and  have  become  an 
intensive  research  area  in  array  signal  processing.  These  techniques  utilize  eigenstructure 
of  the  covariance  matrix.  The  eigenvalues  are  then  used  to  estimate  the  number  of 
sources. 

This  work  considers  that  there  are  D  sources  in  far  field  and  wavefront  is 
considered  planar.  There  are  M  sensors  arranged  in  a  Uniform  Linear  Array.  We  also 
assume  that  number  of  sensors  is  greater  than  number  of  sources.  The  distance  between 
two  sensors  is  d.  Signals  under  consideration  could  be  narrow-band  or  wide-band.  In  the 
narrow-band  case  the  carrier  frequency  is  considered  very  large  as  compared  to  the 
bandwidth  of  the  signal.  If  the  carrier  frequency  is  comparable  to  the  bandwidth  of  the 
signal  then  the  signal  is  considered  as  wide-band  signal.  The  propagation  time  from  one 
sensor  to  another  sensor  is  constant  and  is  generally  approximated  by  a  pure  phase  delay 
in  a  narrow-band  case.  This  approximation  can  not  be  made  whenever  there  is  a  wide¬ 
band  signal.  The  wavefront  impinge  on  the  sensor  with  azimuth  angle#.  An  ULA  is 
shown  in  Figure  1 . 
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We  consider  an  array  of  M  sensors  and  D  point  sources.  Signals  from  these  sources 
impinge  on  the  array  from  directions  ex,e2,...,eD .  The  narrow-band  signals  can  be 
represented  as  time  shifted  version  of  each  other  and  are  expressed  as: 


Xx{t-Tx) 


Y  it)  = 


x2{t-r2) 


XM  ) 


(1.1) 


Where  T|(  is  the  time  it  takes  the  planar  wavefront  to  travel  from  the  source  to  the  origin 
of  the  array  through  the  mth  sensor.  The  signal  is  generally  considered  as  a  complex 
signal  and  can  be  expressed  as 

=  (1.2) 


where  tp  is  an  arbitrary  phase. 
This  could  then  be  expressed  as 


Y(t)  =  x(t-r) 


„-JwT\ 


-iWTi 


„-JWT, 


M 


(1.3) 


The  Mxl  vector  in  the  above  equation  is  generally  referred  as  the  steering  vector.  This 
steering  vector  includes  the  DOA  parameter.  The  signal  model  can  then  be  written  as: 
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'  e-Mi " 

n(t\ 

1 

II 

+ 

n(t)2 

e-imM 

1 

_ 1 

(1.4) 


If  there  are  D  sources  and  M  sensors  then  the  sensor  output  model  can  be  written  as: 


D 

Y (t)  =  £  a(0k)  xk(t-T)  +  n(t) 
k-l 


(1.5) 


where  x/c(t),  i=l,...,D  are  planar  wavefronts  corresponding  to  each  received  source 
signal.  n(t)  is  the  noise  which  is  considered  as  zero  mean  and  white  Gaussian,  a(0)  is  the 
steering  vector 
sin# 


and  T  =  d  - 


(1.6) 


where  d  is  the  distance  between  two  sensors  and  is  considered  as  the  half  of  the 
wavelength. 

Define 


A(tf)=[a(4),  «(»2),  «(«„)] 


a 


te)= 


ax{dk)e 


'  JW0T\  (  @k  , 


au{0k)e 


Jwotm  {@k ) 


(1.7) 


(1.8) 


This  equation  is  also  commonly  known  as  an  array  manifold.  The  output  can  be 
expressed  in  the  matrix  form  as: 

Y(t)  =  A(0)X(t)  +  N(t)  (1.9) 

where 

Y(t),  N(t)  eCM  X(t)  e  C  D  and  A(0)  e  C  MxD 
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Matrix  A  is  called  the  direction  matrix  and  each  of  the  column  vectors  are  the  direction 
vectors  of  the  sources.  Assuming  that  the  noise  is  independent  of  the  signal  with  zero 
mean,  the  covariance  matrix  can  be  expressed  as: 

R  „=E[Y(t)YH(t)]  (1.10) 

If  the  output  array  vector  Y (t)  is  observed  over  K  subintervals  of  duration  AT  seconds 
each,  the  covariance  matrix  can  be  expressed  as  the  snapshot  averaged  cross-product  of 

yk(0  '■ 

d.ii) 

k=i 

1.1  Narrow-band  Multiple  Signal  Classification  algorithm: 

The  Multiple  Signal  Classification  (MUSIC)  algorithm  is  a  high  resolution 
technique.  It  uses  signal  subspace  approach  and  separates  signal  and  noise  subspaces.  It 
finds  the  array  manifold  orthogonal  to  the  noise  space.  Signal  and  noise  subspaces  are 
orthogonal  to  each  other. 

This  means  that  the  eigenvectors  associated  with  the  M  -  D  smallest  eigenvalues  are 
orthogonal  to  the  direction  vectors  making  up  A .  These  observations  form  the  basis  of 
the  MUSIC  algorithm.  We  can  estimate  the  direction  of  arrival  by  finding  direction 
vectors  which  lie  in  the  signal  subspace.  These  vectors  are  the  direction  vectors  that  are 
orthogonal  to  the  noise  subspace.  To  search  the  noise  subspace,  we  form  a  matrix 
containing  the  noise  eigenvectors: 

V„=[qM  -  q„]  (1-12) 
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Since  the  direction  vectors  of  the  incoming  signals  are  orthogonal  to  the  noise  subspace 
eigenvectors,  we  can  say: 

aH  (0)\n\^a(0)  =  O  (1.13) 

where  6  is  the  direction  of  arrival  of  a  signal  component. 

The  direction  of  arrival  can  then  be  estimated  by  finding  the  peaks  of  the  MUSIC 
spectrum  given  by: 


P(0)  = 


a  *  (0)V„VH*a(0) 


-,0e[O,2x] 


(1.14) 


The  D  largest  peaks  in  the  spectrum  will  correspond  to  the  directions  of  arrival  of  the 
signals  impinging  on  the  sensor  array. 

The  MUSIC  algorithm  can  be  summarized  as  follow: 

1.  Collect  input  samples  yk  (t) 

2.  Estimate  the  covariance  matrix  given  by: 

^  t=i 

3.  Compute  eigenvalues  and  eigenvectors  of  the  covariance  matrix. 

4.  Find  number  of  sources  D. 

5.  Find  the  DOA  estimates  by  finding  the  D  largest  peaks  of  the  MUSIC  spectrum 
given  by: 

p{0)=^yvxm~ 
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1.2  Estimating  number  of  sources: 


The  MUSIC  algorithm  depends  on  the  parameter  D,  the  numbers  of  wave  fronts 
impinging  on  the  sensor  array,  for  estimating  the  DOA  of  multiple  targets.  The  Minimum 
Description  Length  (MDL)  was  used  in  order  to  achieve  this  objective  [1-7].  The  MDL  is 
specified  by: 


MDL(Z))=  -2  log 


M 


ru 

k=D+ 1 


( M-D ) 


— S 

-  n  ^ 


M 


M  -  D  *-^k=D+l 

V  J 


+  -(2M  -D)log(N) 


(1.30) 


where  N  is  the  number  of  samples,  k  represents  the  eigenvalues  of  the  covariance  matrix 
and  M  is  the  number  of  sensors.  The  number  of  sources  is  determined  by  finding  a  value 
of  D  which  minimizes  the  MDL  criterion.  The  maximum  number  of  sources  that  can  be 
estimated  is  M  - 1. 


1.3  Wide-band  sources 

Signal  model 

Assume  that  is  a  Uniform  Linear  Array  of  M  sensors  and  there  are  D  wide-band 
sources  with  identical  bandwidth  B  impinging  on  the  array  from  directions ex,e2,...,eD .  In 
the  wide-band  case  the  time  delay  of  planar  wave  propagating  from  one  sensor  to  another 
cannot  be  approximated  as  a  phase  shift.  Assuming  that  the  signals  are  observed  over  a 
finite  interval  T,  we  can  represent  the  signal  xi  by  a  Fourier  series  [1-7]: 

l+m 

xi{t)=TjxiMejWn‘  (L31) 

n=l 

where  Xt  ( wn )  are  the  Fourier  coefficients  given  by: 
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(1.32) 


772 


XiM  =  ^  J  *,(0‘ 


,jw„t 


dt 


-77  2 


2  K 

wn=-^-n,  n=l,...,l  +  m 


(1.33) 


where  w,  is  the  lowest  frequency  and  wl+m  is  the  highest  frequency  included  in  the 

bandwidth  B.  Assuming  that  the  observation  time  T  is  much  greater  than  the  propagation 
delay  across  elements  of  the  array,  we  can  use  a  phase  shift  as  an  approximation  of  the 
time  delay  in  the  Fourier  domain. 

aike~JWnTit  xkM+niM  ( 1  -34) 

k= 1 


The  model  used  to  represent  the  output  vector  is: 

1,(»>*(».)*K)+Nk)  d-35) 

A("'.)  =  [a».(M’.).  a«  (“,  I-  •••  l]  (1-36) 


»*  («’.)  = 


0,(0,  )*-«'*> 


(1.37) 


As  a  result  of  the  Fourier  transform  applied  over  a  time  segment  AT ,  the  array 
output  vector  is  decomposed  into  non-overlapping  narrow-band  components.  The 
covariance  matrix  for  component  wn  can  be  expressed  as: 


Ryy{wn)=E  y{wn)yH{wn) 


This  covariance  matrix  can  also  be  expressed  as: 

Ryy  K  )  =  7!^  K  )  y k  K  ) 

^  fc=l 


(1.38) 


(1.40) 
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Chapter  2 


Review  of  Wide-band  DOA  Algorithms 

There  are  number  of  wide-band  DOA  algorithms  currently  available  in  literature 
mainly  from  IEEE.  Number  of  searches  was  performed  to  locate  wide-band  DOA 
algorithms.  These  searches  were  conducted  at  the  IEEE  digital  library  website  commonly 
known  as  IEEE  Explorer  and  Google.  Searches  at  Google  extend  beyond  IEEE  and  hence 
proved  very  useful.  This  chapter  provides  review  of  wide-band  DOA  algorithms.  Each 
review  contains  introduction,  algorithm  development,  simulation  information  and  our 
conclusions.  Algorithmic  equations  have  modified  as  much  as  possible  to  keep 
uniformity  and  increase  the  readability  of  the  work.  This  work  assumes  that  there  are  D 
wide-band  sources,  a  Uniform  Linear  Array  of  M  sensors,  J  beamformers  (where 
applicable)  and  N  pieces  of  data.  Sources  are  assumed  to  be  in  far  field  so  waves  are 
impinging  on  the  array  as  planar  wave.  Moreover  sources  are  also  Omni-directional. 

Initially  wide-band  DOA  estimation  was  performed  by  estimating  narrow-band 
DOA  at  each  frequency  and  then  combined  to  obtain  the  wide-band  frequency.  This 

approach  is  called  incoherent  approach  [4-5] .  Comparison  of  computational 

requirements  for  all  wide-band  DOA  algorithms  reviewed  in  this  chapter  will  be 
discussed  in  the  next  chapter. 


12 


2.1  Coherent  signal  subspace  method  for  wide-band  sources  by  Wang  &  Kaveh 

It  is  possible  to  combine  the  signal  subspaces  of  different  frequencies  in  order  to 
generate  a  single  subspace  that  will  allow  us  to  determine  the  correct  number  of  sources 
and  directions  of  arrival  [6-7].  The  matrix  R  can  be  used  to  find  the  final  DOA 
estimates.  This  matrix  can  be  formed  by: 

R=STK)R„h)T"K)  <2-‘> 

7=1 

where  J  is  the  number  of  narrow-band  components.  The  matrices  T  are  called 
transformation  matrices  and  can  be  expressed  as 

aiAwo)/aiAwj) 

a2fiM/a2fiM 

aMfAWo)laMp(Wj) 

where  wa  is  the  central  frequency  of  bandwidth  B,  /?  is  the  initial  DOA  value, 
and  ai/}[wj)  is  the  ith  element  of  the  direction  vector  (vw ) . 

The  coherent  signal  subspace  method  for  computation  of  DOA  (wide-band  sources)  as 
proposed  by  Wang  &  Kaveh  [10]  can  be  summarized  as  follow: 

1 .  Collect  data  samples  and  convert  the  samples  into  frequency  domain  using  FFT. 

2.  Estimate  the  covariance  matrix  for  each  frequency  component  given  by: 

Ryy  K ) = 4Z K )  yHk  K ) 

k  k= i 

3.  Compute  eigenvalues  and  eigenvectors  of  the  covariance  matrix. 

4.  Find  initial  estimates  of  direction  of  arrival  by  computing  the  MUSIC  spectrum 
given  by: 
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pip) 


a  ‘MV.V.'M 


5.  Compute  transformation  matrix  focusing  on  central  frequency. 

6.  Compute  eigenvalues  and  eigenvectors  of  the  focus  matrix. 

7.  Find  number  of  sources  D. 

8.  Find  the  final  DOA  estimates  by  finding  the  D  largest  peaks  of  the  MUSIC 
spectrum. 


2.2  Efficient  Wide-band  Source  Localization  Using  Beamforming  Invariance 
Technique  by  Ta-Sung  Lee  (Review) 

A  beamspace  wide-band  source  localization  scheme  is  proposed  which  exploits 
the  concept  of  beamspace  manifold  invariance  [19].  A  principle  of  Least  Squares  (LS)  fit 
is  employed  to  construct  a  beamforming  matrix  for  each  of  the  narrow-band  frequency 
bins  extracted  from  the  wide-band  array  data.  The  beamforming  matrices  perform  the 
same  operation  as  do  the  focusing  matrices  of  Wang  and  Kaveh  [6-7].  In  this  case 
beamforming  matrices  are  chosen  in  such  a  way  as  the  resulting  beamspace  DOA 
matrices  are  essentially  the  same  for  all  frequencies.  The  focused  beamspace  data/noise 
correlation  matrix  pencil  can  then  be  readily  formed  with  the  respective  narrow-band 
beamspace  correlation  matrices  without  any  additional  preliminary  processing. 

A  computationally  efficient  implementation  of  the  beamspace  Root-MUSIC 
algorithm  via  subarray  beamforming  and  banded  transformation  is  developed.  By 
subarray  beamforming,  the  large  order  Root-MUSIC  [1-8]  signal  polynomial  is  first 
reduced  to  one  with  the  order  equal  to  the  beamspace  dimension.  The  algorithm  is 
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further  simplified  by  transforming  the  matrix  representing  the  element  space  noise 
subspace  into  a  banded  form  and  converting  the  reduced-order  signal  polynomial  into 
several  polynomials  with  the  order  equal  to  the  number  of  sources.  These  polynomials  are 
then  rooted  in  parallel  to  determine  the  DOA’s. 

CSSM  Model  Formulation 

This  paper  assumes  that  there  are  D  wide-band  sources  and  M  identical  sensor 
elements  which  have  a  common  passband  width  Bw  centered  at  frequency  fc.  First  of  all 
data  model  needs  to  be  defined  which  uses  following  two  approaches: 

•  The  received  data  is  decomposed  into  J  narrow-band  components  using  a  bank  of 
J  bandpass  filters.  They  are  centered  at  f)  followed  by  the  conventional  I-Q 
demodulation  and  sampling.  The  model  has  J  sets  of  Mxl  complex  array  data 
snapshot  vectors. 

X(rr,f)  =  A(fj  )s(n ;  fj ) +  v(n;  /. )  (2.2) 

where 

s(n;fj)  represents  the  data  received  at  some  reference  point  of  the  array  (Dxl) 
v(n;  fj)  represents  the  noise  present  at  the  M  elements.  (Mxl) 

A(fj)  accounts  for  the  phase  variation  across  the  array  due  to  the  wavefront 
(MxD). 

•  This  model  can  also  be  created  by  forming  narrow-band  components  using 
Fourier  transform  as  done  by  Wang  &  Kaveh  [6]. 

Invariant  means  that  it  is  unaffected  by  the  group  of  mathematical  operations 
under  consideration  and  invariance  means  the  quality  or  state  of  being  invariant. 
Therefore  beamforming-invariance  means  that  beamforming  is  unaffected  by  other 
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factors.  In  a  beam  space  beamforming,  outputs  from  the  array  elements  are  first 
processed  by  a  multiple-beam  beamformer  to  form  a  suite  of  orthogonal  beams.  The 
output  of  each  beam  can  then  be  weighted  and  the  result  combined  to  produce  a  desired 
output.  The  beamformer  can  be  implemented  using  the  FFT.  For  an  M-element  array,  the 
overlapped  orthogonal  beam  can  be  formed  [28]. 

Beamspace  beamforming  requires  a  set  of  beamspace  combiners  to  generate 
weighted  outputs.  The  digital  signal  streams  from  the  antenna  elements  are  passed  to  the 
FFT  processor,  which  generates  K  simultaneous  orthogonal  beams.  The  purpose  of  the 
beam  selection  function  is  to  choose  a  subset  of  these  orthogonal  beams  that  need  to  be 
weighted  to  form  a  desired  output. 

The  approach  in  this  paper  performs  beamspace  transformation  at  each  of  the  J 
frequency  bins  and  choosing  beamforming  matrices  with  some  kind  of  criterion. 
Therefore  beam  patterns  are  identical  for  all  frequencies.  Consider  first  the  patterns 
associated  with  a  single  beam: 

w(F;/.)  =  wffl(F;/.)  (2.3) 

where  j  is  from  1  to  J  and  vv;  is  the  Mxl  complex  weight  vector  employed  at  f  ] . 

It  would  be  difficult  to  find  two  weight  vectors  that  would  produce  identical  beam 
patterns  at  two  different  frequencies.  Due  to  the  discrete  nature  of  the  array,  it  is  not 
possible  to  find  two  weight  vectors  that  produce  completely  identical  beam  patterns  at 
two  different  frequencies.  A  Least  Square  (LS)  fit  method  is  used  as  an  approximation  to 
construct  weight  vectors  that  nearly  provide  frequency-invariant  beam  patterns.  A  simple 
measure  of  the  proximity  between  the  two  patterns  vv(F;/’) and  w(F;/y)  is  the 

generalized  L2  distance: 
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(2.4) 


d.j  =  ^ p(r)\w(r‘,  ft)  -  w(r;  f dr 
a 

dij  =  \p(j)\wf  air-J^-w^air-J^dr  (2.5) 

n 

where  Q  is  a  sector  of  the  unit  sphere  representing  the  Field  of  View  (FOV)  of  the  array. 

The  weighting  function  p(r )  is  enhancing  the  approximation  within  a  pre¬ 
selected  angular  region.  The  weight  vector  computation  for  “beamforming-invariance”  is 
calculated  using  the  following  minimum-distance  criterion 

min|p(F)|wfa(r;/;.)-wfa(F;/y.)|  dr  (2.6) 

Wi.wja 

Above  condition  is  subject  to  that  wt  and  wj  are  not  equal  to  zero.  One 

consideration  that  is  used  as  a  constraint  is  to  obtain  good  SNR  gain  and  low  side  lobes. 
The  weight  selection  process  would  first  require  a  desired  weight  vector  wo  that  would 
be  associated  with  a  pre-selected  reference  frequency  f0  within  the  passband  of  the  array. 
The  Beamforming  Invariance  (BI)  weight  vectors  associated  with  the  J  frequencies  are 
calculated  as: 

ntinj p(r)\w^a(r\fo)- wf a(r;/j)|  dr  (2.7) 

wj  a 

It  may  be  noted  that  f0  need  not  be  one  of  the  frequencies  fj .  Following  solutions  are 
given  by: 

wj=U]1Sjwo  (2.8) 

where 

Uj  =^p(r)a(r,fj)aH(r,fj)dr 
a 
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Sj=\p(r)a{rJj)aH(rJo)dr 

n 

This  technique  is  reducing  the  size  of  data  by  using  a  low-dimension  beamspace 
narrow-band  data  associated  with  /o .  It  is  arguing  that  if  wide-band  source  moves  across 
the  Field  of  View  (FOV)  of  the  array  then  this  set  of  data  would  allow  observation  of 
output  waveforms  at  the  J  beamformer  outputs. 

Beamspace  focusing  with  BI  Transformation 

In  order  to  achieve  effective  reception  of  the  source  signals,  multiple  beams  can 
be  formed  over  the  spatial  band  of  interest.  A  set  of  K  beamforming  weight  vectors  wjk 

are  used  to  simultaneously  form  K  linear  combinations  of  the  array  data  at  fj .  Therefore 

Mxl  element  space  data  snapshot  vectors  are  converted  into  Kxl  beamspace  data 
snapshot  vectors  using  the  following  relation: 

xB(n-fj)=W"x(n-fj)  (2.9) 

where  Wj  are  the  respective  MxK  beamforming  matrices  for  J  frequencies.  K  is  chosen 
such  that  it  is  greater  than  number  of  sources  and  less  than  number  of  sensors.  The 
beamspace  data  snapshot  vectors  will  have  a  similar  form  as  the  original  array  snapshot 
vectors  and  is  given  as: 

xB  («;  fj )  =  B(fj  )s(n;  fj )  +  vB  (n;  f. )  (2. 10) 

where  B(fj)=WJHA(fj)  and  vB(n;fj)=WJHv(n;fj) 

B  and  v  are  the  beamspace  DOA  matrices  and  noise  vectors  respectively. 

The  challenge  would  be  the  selection  of  Wj  and  this  could  be  done  similar  to 
procedure  given  earlier  and  is  also  provided  for  the  convenience. 
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minj P(r)\w?a(r;fJ- wf  a(r;  /,.) |  dr 

wj  n 

It  may  be  noted  that  /  need  not  be  one  of  the  frequencies  /. .  Following  solutions 
are  given  by: 

wj=U]1Sjw0  (2.11) 

where 

U  j  =^p(r)a(r,fj)aH(r,fj)dr 
a 

Sj=^p(r)a(r,fj)aH(r,fJdr 

n 

•  Determine  a  set  of  K  reference  weight  vectors  wok  associated  with  / . 

•  Construct  the  K  weight  vectors  associated  with  the  J  frequencies. 

w]k  =  U:'SjWok 

This  can  be  written  in  matrix  form 

Wj  =u?SjW. 

where  Wo  is  the  reference  beamforming  matrix.  The  beamforming  matrices  for  all  values 
of  j’s  would  have  the  following  form  due  to  the  BI  property: 

B(f  j)  ~  B(f0)  (2.12) 

The  beamspace  data  snapshot  vectors  would  be  fully  characterized  by  a  single 
beamspace  DOA  matrix  B(f0) associated  with/,.  This  will  represent  the  beamspace 

signal  subspace  associated  with  / .  Now  CSSM  approach  can  be  applied  and  obtain  a 
focused  beamspace  data  correlation  matrix. 
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Qxx=Yj  ajEixB  (n;  fj  )x"  (n;  fj ) } 

7=1 

Q,<‘B(QR„Bh(Q  +  Q„ 

7=1 

The  above  operation  is  named  as  BI-CSSM  and  would  produce  an  effective 
narrow-band  beamspace  data/noise  correlation  matrix  pencil  associated  with  associated 
with/,  giving  an  effective  source  correlation  matrix  Rss.  The  following  section  would 
address  the  issue  of  obtaining  beamforming  data  matrix. 

Design  of  Reference  Beamforming  Matrix 

There  are  number  of  publications  and  techniques  for  forming  a  reference 
beamforming  matrix  and  their  references  are  provided  in  this  paper.  A  Chebyshev 
beamformer  can  be  used  which  exhibits  low  side  lobes  for  LES  arrays.  One  of  the  goals 
would  be  to  design  a  beamforming  matrix  that  is  optimal  and  also  has  a  small  focusing 
error. 

Determine  an  orthonormal  basis  for  the  subspace  of  minimum  focusing  error.  Let 
Ew  be  an  MxK’  matrix  satisfying  E".EW  -  I  and  K  <  K  <M  .  Ew  is  used  as  the  reference 

beamforming  matrix  then  from  WJ  - UjlSW0  the  resulting  total  focusing  error  is  given 
by 

1  /  2 

Ef  II  dr=tr{E*SvEw)  (2.13) 

J  j= i 

and 

Sv=\i  (S0-S“U:'Sj) 

J  7=1 
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Qxx=Yj  ajEixB  (n;  fj  )x"  (n;  fj ) } 

7=1 

Q,<‘B(QR„Bh(Q  +  Q„ 

7=1 

The  above  operation  is  named  as  BI-CSSM  and  would  produce  an  effective 
narrow-band  beamspace  data/noise  correlation  matrix  pencil  associated  with  associated 
with/,  giving  an  effective  source  correlation  matrix  Rss.  The  following  section  would 
address  the  issue  of  obtaining  beamforming  data  matrix. 

Design  of  Reference  Beamforming  Matrix 

There  are  number  of  publications  and  techniques  for  forming  a  reference 
beamforming  matrix  and  their  references  are  provided  in  this  paper.  A  Chebyshev 
beamformer  can  be  used  which  exhibits  low  side  lobes  for  LES  arrays.  One  of  the  goals 
would  be  to  design  a  beamforming  matrix  that  is  optimal  and  also  has  a  small  focusing 
error. 

Determine  an  orthonormal  basis  for  the  subspace  of  minimum  focusing  error.  Let 
Ew  be  an  MxK’  matrix  satisfying  E".EW  -  I  and  K  <  K  <M  .  Ew  is  used  as  the  reference 

beamforming  matrix  then  from  WJ  - UjlSW0  the  resulting  total  focusing  error  is  given 
by 

1  /  2 

Ef  II  dr=tr{E*SvEw)  (2.13) 

J  j= i 

and 

Sv=\i  (S0-S“U:'Sj) 

J  7=1 
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Design  Example 

An  array  manifold  vector  associated  with  an  LES  array  consisting  of  M  identical 
elements  and  operating  with  frequency/is  given  by: 

a{u\f)  =  [l,em,ej2fu, . (2.14) 

where  t  -  2nd  / c  and  d  is  the  spacing  between  two  adjacent  elements.  The  sine-space 
angle  u  is  defined  as  u=sin(0)  and  60s  the  angle  measured  with  respect  to  the  broadside 
of  the  array.  First  sensor  element  is  assumed  to  be  the  reference  point  of  the  array. 

The  frequency  band  is  decomposed  into  J=33  uniformly  distributed  subbands.  The 
reference  beamforming  matrix  Wo  is  constructed  at  /  =  f\ .  by  first  formulating  Wd  and 
using  following  equations: 

min  IK -W',,  f 

Wo=EwV 

WQ=Ew{E:EwYlEHwWd 

Use  K’=9  and  Wd  is  composed  of  weight  vectors  associated  with  k=7  Chebyshev 
beams  with  -30  dB  sidelobles  point  at  0,  ±  7.7,  ±  15.5  and  ±  23.6  in  degrees.  To  alleviate 
the  grating  lobe  problem  fo  =  /j  is  chosen.  Remaining  thirty  two  beamforming  matrices 

are  computed  via  W-  =U~j1SjW0,.  The  FOV  was  set  to  be  -1.0  to  1.0  in  u  domain.  The 

weighting  function  was  chosen  as  1.0  and  0.5.  These  assumptions  were  verified  by 
computing  the  normalized  focusing  error  spectrum. 

BI-CSSM/ROOT  -Form  Eigen-Based  DOA  Estimation 

Now  consider  only  the  operation  associated  with  fo  and  omit  the  argument  of 
frequency  in  the  relevant  terms.  The  beamspace  eigen-based  methods  can  be  applied  on 
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the  BI_CSSM  focused  data  dictates.  It  is  known  that  the  source  DOA’s  should  be 
determined  via  the  null  spectrum 

<D(F)  =  aH(r)W0EBPEBHW0Ha(r )  (2.15) 

where 

Eb  is  referred  to  as  the  beamspace  noise  eigenvector  matrix  which  is  Kx(K-D).  It 
contains  generalized  eigenvectors. 

P  is  a  positive  semi-definite  matrix  serving  to  weight  the  respective  columns  of  EB . 

The  null  spectrum  is  converted  into  the  2(M-l)th  order  signal  polynomial 
d>(z)  =  aT(Z-lW0EBPE”W0Ha(z) 
where  a(z)=[l,z,. .  ..zM_1]T  and  z  -  ej^°u . 

There  are  D  signal  roots  z\  extracted  from  d>(z)  .  Their  DOA’s  can  be 
determined  by  ut  =  arg{f(  where  i  is  from  1  to  D.  These  coefficients  of  d>(z) 

exhibit  conjugate  symmetry  such  that  the  corresponding  2(M-1)  roots  from  M-l 
conjugate  reciprocal  pairs.  As  a  result  only  M-l  distinct  values  are  observed  regarding  the 
phase  angles  of  2(M-1)  roots. 

The  root-form  methods  may  be  computationally  expensive  due  to  the  need  of  large 
order  polynomial  rooting.  A  new  method  is  proposed  which  may  be  more 
computationally  efficient  and  has  following  three  stages: 

1.  The  signal  polynomial  is  reduced  from  order  2(M-1)  to  2(K-1)  via  judiciously 
performed  subarray  beamforming. 
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2.  The  reduced  order  signal  polynomial  is  converted  into  several  2Dth  or  (Dth)  order 
polynomials  via  a  banded  transformation  of  the  corresponding  reduced  noise  EV 
matrix. 

3.  Signal  roots  are  extracted  by  rooting  these  polynomials  in  parallel. 
Polynomial-Order  Reduction  via  Subarray  Beamforming 

Consider  an  M-element  LES  array  which  consists  of  L=M-K+1  and  overlaps  K- 
element  subarrays.  Define  the  KxM  selection  matrices  that  select  from  the  full  array  data 
snapshot  vectors  the  respective  subarray  data  snapshot  vectors 

T ]x(n)  =  x(^\n) 

where  x(^(n)  denotes  the  Kxl  data  snapshot  vectors  received  at  the  1th  subarray  and  T)  is 
given  by 

=  [0KX(l- 1)  I  I KXK  I  OkX(M-K-1+ 1)] 

Subscripts  indicate  the  sizes  of  the  respective  identity  and  zero  matrices.  Assume  that  the 
same  Kxl  weight  vector  g  is  applied  at  each  of  the  subarrays  that  would  produce  a  set  of 
Lxl  vectors. 

V4V)~ 

gHxf(.n ) 

xL(n)=  .  =GHx(n )  (2.16) 

_gH4L\n)_ 

It  represents  the  data  snapshots  received  at  the  L  subarray  beamformer  outputs. 
xL(n)  is  the  data  snapshot  vector  of  L-element  LES  array  which  is  also  denoted  as  AL  and 
apply  an  Lxl  vector  c  to  form: 

xB(ri)  -cH xL{n)  —  (Gc)H x{ri)  —  wH x(n)  (2.17) 
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It  can  be  seen  that  w=Gc  and  is  the  effective  weight  vector  acting  on  the  entire 
array.  We  may  consider  the  beamforming  to  be  performed  first  on  Al  and  treating  the  L 
subarray  as  super  elements.  We  may  also  apply  an  Lxl  weight  vector  c  to  produce  a  set 
of  Kxl  vector  data  snapshots. 

xK(n)  =  Sc'*xx(n)  =  CHx(n) 

1=1 

L 

where  C  -  ^ 

/= i 

These  vector  data  snapshots  can  be  regarded  as  “data  snapshot  vectors”  obtained 
from  Ak-  Applying  a  Kxl  weight  vector  g  forms: 

xB(n)  =  gH  xL(n)  =  ( Cg)H  x(n )  =  wH  x(n)  (2.18) 

It  can  be  concluded  that  the  full  weight  vector  exhibits  two  types  of  decomposition: 

W=Gc=Cg 

Using  the  banded,  Toeplitz  structure  of  C  or  G 

wHa(z)  =  {gHak(z)}{cHaL(z)} 

Similarly  we  can  get 

aT(z~l)w  =  {aTK(z~l)g}{aTL(z~l)c} 

A  beamspace  transformation  scheme  based  on  the  concept  of  weight  vector 
decomposition  can  be  developed.  K  reference  weight  vector  can  be  constructed  using  the 
following  form: 

W0 k  ~  GkC  ~  Cgk 

The  w  contains  the  same  weight  factor  c.  Now  putting  in  matrix  form  and  using  above 
equations: 
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WHa(z)  =  {GHak(z)}{cHaL(z)} 
aT(z-l)W0={aTK(z-l)G}{al(z-l)c} 


Substituting  these  equations  in  <t>(z)  =  aT(z  1  )W0EBPEBW^  a(z)  yields 

^{z)^{al{z-l)ccHaL{z)}{al{z-x)GEBPE^GHaK{z)}  (2.19) 

It  can  be  seen  that  this  is  decomposed  into  two  individual  factors  accounting  for  c 
and  G.  The  factor  involving  c  is  known  a  priori  and  thus  there  is  no  information  about  the 
DO  As.  The  DOA  estimates  can  be  determined  with  the  2(K-l)th  order  polynomial. 

(Z)  =  4  (z"‘  )GEBPE»GHaK  (z) } 

This  equation  suggests  that  this  is  the  signal  polynomial  associated  with  K- 
element  LES  array  generated  by  the  noise  EV  matrix  GEB .  This  equation  leads  to 
substantial  reduction  in  computational  load  if  K«M. 

Parallel  Processing  via  Banded  Transformation 

Subarray  beamforming  somewhat  simplifies  the  polynomial  rooting.  The 
difficulty  remain  with  the  procedure  of  choosing  D  signal  roots  out  of  the  2(K-1)  roots  of 
<i> K(z)  ■  A  scheme  has  been  proposed  by  others  to  convert  the  DOA  estimation  problem 
into  that  of  rooting  a  Dth  order  polynomial.  They  exploit  the  fact  that  the  ideal  noise 
subspace  associated  with  an  LES  array  is  spanned  by  the  columns  of  a  banded  Toeplitz 
matrix  with  bandwidth  D+l. 


(2.20) 


Performing  some  algebraic  manipulation  which  leads  to: 
®K(z)  =  aTD+l(z-l)FD(z-1)TBPTBHD(z)FHaD+l(z)} 
where  a  D+1(z)=a(z)  and  D(z)  is  the  diag  { l,z, . ,zK  D+1 } 
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The  polynomial  is  associated  with  a  (D+l)  element  LES  array  except  for  the  z-dependent 
term  D  (z).  In  order  to  fully  exploit  this  reduction  in  dimension,  replace  D  (z)  with  a 
constant  matrix  by  fixing  z=z0. 

I  Zo)  =  aTD+1(z-1)FD(z0-l)TBPTBHD(z0)FHaD+l(z)  (2.21) 

Comparing  two  <5? K(z)  equations,  it  can  be  seen  that  in  order  to  achieve  the 
performance  of  working  with  original  signal  polynomial,  we  must  choose  z0  ~  z,  =  ejf°Ui 

in  estimating  Ui.  It  indicates  that  a  set  of  reduced  order  signal  polynomial  be  constructed 
with  Zj,  i=l  D.  The  problem  with  this  approach  is  that  it  requires  the  knowledge  of 
the  DOA’s  that  is  being  estimated.  The  effective  spatial  passband  of  the  reference 
beamforming  matrix  can  be  decomposed  into  Is  disjointing  sectors  centered  at  um, 
m=l,. . ..,  Is  and  construct  the  following: 

<Mz  I  ZB)  =  al+l(zAFD(zm-l)TBPTBHD(zm)FHaD+l(z)  (2.22) 

with  zm  -  and  m=l,....,  Is .  Q  K(z)  can  be  approximated  with  the  above  equation 

for  the  mth  sector  with  high  accuracy  for  a  small  sector  size.  This  equation  can  be  rooted 
in  parallel  and  obtain  Is  set  of  roots  of  zim. 

Roots  need  to  be  picked  and  we  should  pick  roots  which  are  closer  to  the  unit 
circle.  This  DOA  estimator  is  suboptimum  according  to  the  author.  It  would  also  result  in 
degradation  in  estimation  accuracy  at  low  SNR.  Even  with  the  true  DOA,  the  DOA 
estimates  may  not  be  exactly  identical  as  signal  roots  may  not  lay  on  the  unit  circle.  The 
parallelized  estimator  behaves  as  a  mixture  of  the  root-form  and  spectral-form  estimator. 

We  note  that  under  no  noise/error  conditions  is  the  banded  matrix  is  Toelplitz.  We 
may  assume  that  F  is  approximately  rank  one.  This  could  be  done  using  following 
replacement: 
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(2.23) 


FD{zm-l)TBPT» D(zJFh aD+l(z)  =>  FD(zm-l)TBPmP»TBH D(zm)FH 

where  Pm  is  a  (K-D)xl  vector.  Subscript  m  is  used  to  emphasize  the  dependence  on 
different  sectors.  With  the  above  structure  we  need  only  work  with  the  set  of  Dth  order 
polynomials: 

<S>Az\zJ  =  PmHTBHD(zm)FHaD+l(z) 

Algorithm  Summary 

•  Construct  beamforming  matrices  and  the  BI  transformations. 

•  Perform  a  KxK  generalized  eigen-decomposition 

•  Solve  in  parallel  K-D  systems  of  equations  of  size  K-D. 

•  Rooting  in  parallel  Is  2Dth  (or  Dth)  order  polynomials. 

Design  of  Subarray-Based  Reference  beamforming  Matrix 

The  factorization  of  the  equation  w0k  -  Gkc  -  Cgk  does  not  hold  for  a  particular 

desired  beamforming  matrix  W)/.  One  way  to  retain  the  merits  of  Wj  using  subarray 
beamforming  is  to  selectively  choose  c  and  gk  so  that  Wo  is  close  to  W(j.  Using  the  LS  fit 
technique  leads  to  the  following  problem: 

min  min 

r - N  2  t - A~- - v  K '  ^ 

c,gl,...,gkfVo-Wd  \F=c,gx,..gkY}w„k-wdk  |  (2.24) 

k= 1 

Invoking  the  structure  of  w0k  -  Gkc  -  Cgk  we  can  re-write  above  equation  as: 


C,G  CG-Wd 


=  C,G 


Gc-w * 


|2 


This  equation  has  no  closed  form  solution  in  general.  The  problem  should  be 
decomposed  into  two  individual  stages  for  which  in  one  stage  we  solve  for  the  common 
weight  factor  c  whereas  in  the  other  we  solve  for  the  uncommon  weight  factor  gk. 
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Assuming  that  the  initial  guess  of  c  is  available  and  solving  the  left  hand  side  of  G  we 


get: 


d  =  (cHcylcHwd 


Constructing  G  and  then  solving  the  right  hand  side  for  c  yields: 

c  =  (GHG)~lGHwd 

Now  construct  C  and  solve  for  new  G.  This  procedure  is  then  alternatively  executed  until 
the  solution  converges. 

Summary 

The  algorithm  decomposes  data  using  bandpass  filters  into  J  frequency  beams.  It 
performs  beamspace  transformation  and  computes  weights  using  least  square  method.  It 
then  computes  beamspace  data  matrix  and  focuses  on  single  reference  frequency  which 
would  be  something  similar  to  CSSM  method.  It  performs  transformation  into  K 
beamspaces  and  forms  beamspace  data  matrix.  This  beamspace  data  matrix  then  focused 
on  a  single  reference  frequency  out  of  J  frequency  bands.  The  design  of  beamspace  data 
matrix  is  described  which  requires  first  design  of  beamforming  matrices.  They  are  again 
designed  in  a  least  square  sense.  The  problem  is  then  reduced  to  beamspace  data 
correlation  and  noise  matrices.  Authors  then  apply  their  own  derived  root  MUSIC 
algorithm  which  could  be  substituted  with  the  MUSIC  algorithm.  The  algorithm  does  not 
require  any  preliminary  DOA  estimates. 

This  DOA  estimator  is  suboptimum  according  to  the  author.  It  would  also  result 
in  degradation  in  estimation  accuracy  at  low  SNR.  Even  with  the  true  DOA,  the  DOA 
estimates  may  not  be  exactly  identical  as  signal  roots  may  not  lay  on  the  unit  circle.  The 
proposed  parallelized  estimator  behaves  as  a  mixture  of  the  root-form  and  spectral-form 


29 


estimator.  This  approach  provides  an  excellent  alternative  to  CSSM  but  frequency 
decomposition  and  calculations  of  weight  and  beamspace  data  matrices.  There  may  be 
some  alternative  ways  to  calculate  weight  and  compute  data  beamspace  matrices.  In  the 
end  it  again  uses  something  similar  to  MUSIC  algorithm  to  compute  DOA. 

2.3  Multiple  Broad-Band  Source  Location  Using  Steered  Covariance  Matrices 
Jefrey  Krolik  and  David  Swingler  (Review)  [17] 

This  paper  is  based  on  space  time  statistics  called  the  Steered  Covariance  Matrix 
(STCM).  It  is  obtained  by  measuring  the  covariance  of  the  time-domain  array  outputs 
after  delays  have  been  inserted  to  steer  a  conventional  Delay-and-Sum  (DS)  beamformer 
beam  [17].  This  technique  avoids  source  localization  problem  by  defining  a  broad  band 
covariance  matrix  having  a  rank  one  characterization  regardless  of  source  spectral 
content  or  source  location.  The  drawback  of  this  approach  is  that  it  is  compute  intensive. 
Model  Formulation  and  the  steered  covariance  matrix 

We  consider  an  array  of  M  wide-band  sensors  and  D  wide-band  point  sources. 
The  output  of  a  sensor  located  at  xm  is  denoted  by  ym(t)  observed  over  a  time  interval  of 
T  seconds  and  is  expressed  as: 

y.m  =  i,u,(t-T„m>  +  vJO  (2.25) 

(=1 

where  Uj(t ),  i=l,...,D  are  stationary,  zero-mean,  random  processes  corresponding  to  each 
received  source  signal  and  vm(t)  is  the  noise.  zm  (0! )  is  the  signal  propagation  delay  to 

the  mth  sensor  when  a  source  is  located  at  8j . 
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The  locations  0i  are  the  parameters  that  need  to  be  estimated  from  a  finite-time 

observation  of  the  sensor  outputs.  The  STCM  can  be  defined  using  the  complex  analytic 
representation  of  the  sensor  outputs.  The  delay  and  sum  beamformer  enhances  the 
reception  of  signals  emanating  from  location  0  by  inserting  a  delay  of  Tm  (0)  at  the 
output  of  each  sensor.  The  complex  analytic  representation  of  the  delay  and  sum 
beamformer  output  b  (t,  0)  with  steering  in  direction  <9  is  formed  by  taking  the  weighted 
sum  of  the  elements  of  y  (t,  0)  and  is  expressed  as: 

b(t,0)  =  wTy(t,0) 

where  w’s  are  a  vector  of  real- valued  array  shading  weights.  The  expected  beam  power 
would  be: 

z(0)  =  E{\b(t,0)\2} 

This  could  then  be  expressed  as  z( 0)=  wT R( 0 )w 

where  R  (0)  =E {y(t,  0)  y(t,  0)H)  is  the  Steered  Covariance  Matrix  (STCM)  and  is 
corresponding  to  direction  0 . 

The  well  known  Cross  Spectral  Density  Matrix  (CSDM)  is  a  function  of  temporal 
frequency  and  its  corresponding  weighting  vector  depends  both  on  the  array  shading  and 
steering  angle.  R  ( 0  )  is  a  function  of  steering  0  and  w  is  a  constant  shading  vector.  The 
main  advantage  of  using  R  ( 0  )  in  wide-band  case  is  that  the  number  of  statistical  degrees 
of  freedom  available  to  estimate  R  ( 0  )  is  approximately  equal  to  the  time-bandwidth 
product  of  the  sources  rather  than  the  usually  much  smaller  number  of  snapshots  used  in 
estimating  each  narrow-band  CSDM. 
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Structure  of  R  (0  )  can  be  used  for  a  wide  variety  of  high  resolution  source 
location  methods.  Consider  a  Uniform  Linear  Array  (ULA)  with  M  sensors  at  d  distance 
apart  at  positions  xm=md,  m=0,. . .  ,  M-l. 

Tm(0)  =  rn.T(0)  (2.26) 

where  t(0)  -  d  /  csin(#) ,  c  is  the  propagation  speed  and  0  is  the  bearing  of  the  source 

D 

relative  to  array  broadside.  Using  the  model  of  ym(t)-'^_iui{t-Tm{0i))  +  vm{t),  the  jkth 

(=1 

element  of  R  ( 0  )  is  a  function  of  m=j-k  given  by: 

rjk  (0)  =  2  Pi  (m.(T(0)  -  m ))  +  t]m  ( m.T(0 ))  (2.27) 

i—X 

where  pi  (t)  is  autocorrelation  function  of  the  ith  source  and  ijm  (t)  is  the  cross-correlation 

function  of  the  noise  received  at  the  array  coordinate  origin  and  the  mth  sensor. 

When  the  steering  direction  is  the  direction  of  the  source  and  the  source  is  aligned 
with  the  steering  direction  then  the  STCM  contains  a  constant.  This  perfectly  coherent 
component  is  equal  to  the  source  power  regardless  of  its  spectral  signature.  One  approach 
to  determine  the  power  level  of  a  point  source  in  the  steering  direction  0  is  to  accurately 
estimate  the  level  of  the  dc  constant  term.  This  estimation  is  performed  for  a  closely 
spaced  set  of  steering  directions  yields  a  wide-band  spatial  power  spectral  estimate.  It  can 
also  be  noted  that  DS  beamforming  corresponds  to  making  a  conventional  Blackman- 
Tukey  estimate  of  the  dc  component  of  r  (m,  0).  Later  on  minimum  variance  and  linear 
predictive  spectral  estimation  are  examined  as  methods  of  estimating  the  dc  component 
for  each  steering  direction. 
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Estimation  of  the  Steered  Covariance  Matrix  (STCM) 

An  estimate  of  the  Steered  Covariance  Matrix  (STCM)  is  obtained  using  a  simple 
relationship  between  the  STCM  and  CSDM.  Statistics  of  the  estimated  STCM  can  be 
expressed  in  terms  of  the  Wishart  characteristic  function.  Complex  time  domain  vector  of 
M  sensor  outputs  can  be  expressed  as  [y^  (t),  ycx  ( t ycM_x  ( t ) }r  over  the  time  interval  (- 
T/2,T/2)  in  terms  of  the  frequency-domain  vectors  as 

Y(k)  =  [Y0(k),Yx(k),...YM_x(k)]T  (2.28) 

With  elements  Ym(k)  corresponds  to  the  Fourier  series  coefficients  of  ym(t)  at  frequency 
6)k  =  2 7ik  /  T  The  sensor  outputs  are  approximately  band-limited  to  0),  <  co  <  coh .  The 
steered  sensor  output  vector  y(t,  6)  can  be  expressed  as: 

y(t,0)  =  fjTk(0)Y(k)ej^t  (2.29) 

k=l 

where 


ejCOkTo(0 ) 

Q  ej0Jkrl{6) 


Tk{d)  = 


(2.30) 


0 


gjGkTM-\(9) 


Substituting  above  equation  of  y(t,  0)  into  R  (0)  and  under  the  assumption  of  large  T, 
the  STCM  can  be  expressed  as: 

R(0)  =  (d)K(cok  )Tk  0 9)H  (2.31) 

k=l 

where  K(cok )  =  E{Y(k)Y(k)H }  is  the  conventional  un-steered  CSDM  at  frequency  cok . 
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This  R  (0)  matrix  is  similar  to  coherently  focused  covariance  matrix  of  CSS 
method  proposed  by  Wang  and  Kaveh  [2]  for  the  case  where  all  sources  in  the  field  are  in 
a  single  group  unresolved  by  a  conventional  delay  and  sum  beamformer.  In  STCM  the  R 
( 6 )  matrix  is  computed  for  each  steering  angle  6  making  it  computationally  intensive. 
However,  it  eliminates  the  source  location  bias  resulting  from  errors  made  in  forming 
focusing  matrices  required  by  the  CSS  method. 

The  relationship  between  matrices  K  and  R  {6)  as  given  in 

h 

R(6)  -  jTk  (6)K(cok  )Tk  (0) H  suggests  a  natural  way  of  estimating  R  ( 0 )  by  using  finite- 

k=l 

time  CSDM  estimates  for  K .  A  common  method  of  forming  K  from  discrete-time  sensor 
outputs  is  to  divide  the  T  second  observation  into  N  non-overlapping  segments  of  AT 
seconds  each  and  then  apply  the  FFT  to  obtain  uncorrelated  frequency  domain  vectors 
Yn( k)  for  each  segment  n=l,...,N.  The  cross  spectral  density  matrix  at  each  frequency 
cok  is  then  estimated  and  then  substituted  to  estimate  of  R  ( 0 )  matrix.  An  attractive 

feature  of  the  STCM  estimate  given  is  that  its  statistics  can  be  expressed  in  terms  of  the 
Wishart  characteristic  function. 

Steered  Covariance  Source  Location  Methods 

It  is  derived  by  finding  the  beamformer  weight  vector  w  which  minimizes  the 
beam  power  given  by  z(  0 )=  wT R(  0 )w  subject  to  the  constraint  that  the  processor  gain  is 
unity  for  a  broad-band  plane  wave  in  direction  0 .  The  problem  is  viewed  as  one  of 
estimating  the  dc  component  of  the  STCM  steered  in  direction  0  by  means  of  a 
Minimum  Variance  (MV)  approach.  In  either  case,  this  technique  has  the  effect  of 
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choosing  w  to  minimize  the  power  contribution  from  sources  and  noise  not  propagating 
from  direction  0 .  The  resulting  STCM-based  spatial  spectral  estimate  is  given  by 
Zsmv  =  |l"  K(<9)"'ll  1  (2.32) 

where  1  is  an  Mxl  vector  of  ones. 

A  finite  estimate  of  ZStmv  can  be  obtained  by  substituting  the  estimate  of  R  (0) 

h 

matrix  similarly  estimated  by  R(0)  =  ^Tk(0)K(a>k)Tk(0)H  .  A  complete  broad-band 

k=l 

spatial  power  spectral  estimate  is  then  formed  by  computing  Zsmv(0 )  for  a  set  of 
steering  angles  which  span  the  locations  of  interest. 

For  wide-band  sources,  the  Minimum  Variance  Distortionless  Response  (MVDR)  beam 
power  is  obtained  by  summing  narrow-band  beam  powers  over  the  band  of  interest  and  is 
given  as: 

z„„*  (»)  =  Z[  A  (9  f  Kia,  )-■  D,  (0)]  (2.33) 

k=l 

The  steps  in  the  STMV  method  are  as  follows: 

1 .  Form  estimated  Cross-Spectral  Density  Matrices  k  over  the  frequency  band  of  interest. 

K(d>l)  =  ^Y.(k)Y„(k)“  (2.34) 

n= 1 

2.  Compute  estimated  steered  covariance  matrices  for  each  steering  direction  0  of 
interest. 

R{0)  =  (0)K(cok  )Tk  (0)H  (2.35) 

k-l 

3.  Compute  R  (0)  in  the  following  equation  and  form  Z  for  each  steering  direction  0  to 
obtain  a  broad-band  spatial  power  spectral  estimate. 
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(2.36) 


Z 


STMV 


[lw  R(6)~l  1]_1 


Compute  linear  prediction  coefficients  from  the  steered  sensor  outputs  and 
approximate  or  model  them  for  the  Autoregressive  (AR)  model.  Compute  forward  and 
backward  prediction  error  sequence  and  linear  prediction  coefficients.  The  linear 
prediction  coefficients  which  minimize  forward  and  backward  prediction  error  sequence 
are  simply  complex  conjugates  of  each  other.  Exploiting  this  property  the  Steered  Linear 
Prediction  (STLP)  method  is  simply  the  minimizing  the  sum  of  forward  and  backward 
squared  prediction  errors.  The  steered  linear  predictive  estimate  of  the  spatial  power 
spectrum  in  direction  6  is  given  by 


PM 

i -ix.(0> 


(2.37) 


This  is  also  known  as  a  maximum  entropy  spectral  estimate  evaluated  at  dc.  The 
computation  of  above  equation  and  as  the  minimization  of  average  forward  and  backward 
prediction  squared  error  can  be  obtained  as  the  solution  to  the  following  equation. 


RLmL(0) 


pL(a) 

0L 


Summary  of  the  STLP  algorithm 

•  Form  spatially  averaged  cross  spectral  density  matrix  estimates  K  over  the 
frequency  band  of  interest. 

•  Compute  spatially  averaged  steered  covariance  matrices  R  for  each  steering 
direction  6  of  interest. 
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to  obtain  vector 


Solve  the  augmented  normal  equation  RL(0)aL(0)  = 


pL(a) 

0L 


STLP  coefficients  aL(0)  and  p,{9)  for  each  steering  direction  9 . 


Compute  equation  Zstl  (0)  - 


£M 

l-Z 

m-\ 


Simulation  Study 

This  study  simulates  using  a  ULA  with  16  sensors  with  inter-element  spacing  of 
7.5  m  corresponding  to  a  half  wavelength  frequency  of  100  Hz  is  considered.  The 
location  parameter  6  is  the  bearing  of  the  source  relative  to  endfire.  The  Rayleigh  limit 
of  angular  resolution  for  this  array  is  approximately  2/(M-l)=0.133  radians  or  7.62 
degrees.  Two  source  signals  were  modeled  as  temporally  stationary  and  mutually 
uncorrelated  zero  mean  Gaussian  processes  with  bandpass  auto-spectra. 

Si(f)  Li  forfo-BW/2  <f<f0+BW/2 

0  otherwise 

where  each  source  has  the  same  frequency  fo=100Hz  and  bandwidth,  BW=40Hz.  The 
noises  vm(t),  m=l..M  were  taken  to  be  stationary  and  mutually  uncorrelated  zero-mean 
Gaussian  bandpass  processes  independent  from  the  signals  and  watch  with  the  auto¬ 
spectrum: 

Sv(f)  Lv  forfo-BW/2  <f<f0+BW/2 

1  otherwise 

defined  over  the  same  frequency  band  as  the  source  signals.  The  SNR  of  the  ith  source 
denoted  SNRj  is  defined  from  the  above  by  SNRi=Li/Lv  Estimated  cross-spectral  density 
matrices  were  formed  using 
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(2.38) 


«®*>  =  T7  EWW" 

^  ’  n=l 

From  sensor  outputs  where  a  segment  length  AT  =0.8  seconds  was  used.  For  each 
segment  the  array  output  is  decomposed  into  B=h-1+1=33  narrow-band  DFT  bins.  The 
number  of  segments  or  snapshots  N  employed  to  estimate  each  k  is  varied  in  this 
simulation  study. 

Summary 

We  have  used  our  simulation  program  with  previously  generated  data  outlined  in  [6]. 
Our  program  did  not  give  two  peaks  as  expected.  Problems  could  be  in  two  different 
areas: 

1.  There  may  be  some  missing  information  regarding  transformation  matrix  that 
could  lead  to  not  getting  appropriate  result. 

2.  The  paper  uses  random  processes  as  data  which  may  contribute  to  getting  two 
peaks  in  their  case  and  not  in  our  case. 

3.  The  third  conclusion  we  have  is  that  this  technique  requires  computation  for  each 
steering  direction  0  increasing  the  computation  requirements  but  eliminates  the 
selection  of  initial  focusing  angle. 

2.4  Focused  Wide-Band  Array  Processing  by  Spatial  Re-sampling  by  Jefrey  Krolik 
and  David  Swingler  (Review)  [18] 

This  focusing  technique  reduces  each  wide-band  source  in  multigroup  scenarios 
to  a  rank  one  representation.  This  approach  does  not  require  preliminary  estimates  of  the 
source  locations.  The  method  is  based  on  simply  adjusting  the  spatial  sampling  rate  or 
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“spatially  re-sampling”  the  array  outputs  as  a  function  of  temporal  frequency  so  that 
wide-band  sources  are  aligned  in  the  spatial  frequency  domain.  [18] 

In  this  case  an  output  of  a  discrete  array  of  M  sensors  is  considered  as  the  result  of 
spatially  sampled  a  continuous  linear  array.  Let  y(x,  co)  denote  the  field  incident  upon  a 
line  array  positioned  along  the  x  axis  at  temporal  frequency  co.  The  y(x,  co)  is  expressed 
as: 

y(x,  ry)  =  X  S,  (co)e’C0a-x  +  v(x,  co)  (2.39) 

(=1 

where  S  is  the  temporal  Fourier  transform 
ai  -  sin(#( )  /  c  is  the  slowness 
6  is  the  bearing  angle. 

Assume  a  Uniform  Linear  Array  (ULA)  with  M  sensors  is  spaced  at  a  uniform 
distance  d  meters  apart.  The  yd  (m,  co)  is  the  output  of  the  sensor  located  at  x=md  and  is 
expressed  as: 

D 

yd  (x,  cq)  =  Yj  si  (<n)e]Coa'md  +  v(md,co )  (2.40) 

i= 1 

A  narrow-band  covariance  matrix  R(co)  can  be  computed  from  sensor  array  outputs.  This 
narrow-band  covariance  matrix  R  (co )  can  be  expressed  as: 

R{co)  -  A{co,a)Ps{co)A{co,a)H  +  Rv(co)  (2.41) 

where  A  is  the  MxD  source  direction  matrix  of  the  ith  source,  P  is  an  unknown  DxD 
source  spectral  density  matrix  Rv  is  the  noise  matrix. 

The  goal  of  focusing  methods  is  to  transform  the  sensor  outputs  yd  (m,  co)  in  such  a  way 
that  it  results  in  a  source  direction  matrix  which  is  constant  for  all  frequencies  within  the 
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common  bandwidth  of  the  sources.  A  direction  matrix  A  is  obtained  for  all  co  by 
adjusting  the  spatial  sampling  interval  d  as  a  function  of  temporal  frequency  co.  A 
focusing  frequency  coo  is  also  selected. 

The  frequency  dependent  spatial  sampling  interval  is  denoted  by  d  (co).  The 
wide-band  focusing  via  spatial  re-sampling  can  be  achieved  by  letting  d  (ta)  =d  coj  t a. 
If  we  substitute  d  ( ta )  in  the  above  expression  we  get  following  focused  array  outputs: 

yd(x,t0)  =  fjS,(toy’"'~'  +v(^t,®)  (2.42) 

,=i  0) 

In  order  to  avoid  any  spatial  aliasing,  d  (co)  must  be  chosen  such 
that  <JJ/  c  <  7l /  d(co) .  This  implies  6)0<7ie/d  relation  should  be  true.  It  is  known  that  the 

sensor  spacing  should  be  half  wavelength  of  the  highest  source  frequency  tamax  then  this 
will  result  in  focusing  frequency  0)0  should  be  less  than  highest  source  frequency  giving 
oj0  <  tamax .  The  re-sampling  of  y(x,  co )  from  the  outputs  of  a  discrete  line  array  yd(m,  co ) 
at  a  finite  set  of  frequencies  con  (  where  n  is  from  1  to  B)  should  ensure  co{)  /  con  =  Ln/  Kn. 
Therefore  the  sampling  interval  should  follow  this  relation 

d(  con )  =d.Ln/Kn 

Therefore  spatially  re-sampling  the  data  at  frequency  con  involves  changing  the 
spatial  sampling  rate  by  a  factor  of  Kn/Ln.  This  is  also  achieved  by  selecting  the  focusing 
frequency  co0  equal  to  the  minimum  frequency  of  the  sources  CDmm  .  If  elements  are  spaced 

a  half  wavelength  apart  at  <wmax  ( d  =  tvc  1 0)ma  )  this  implies  Ln  <  Kn  for  all  n.  Hence  re¬ 
sampling  becomes  an  interpolation  by  a  factor  of  Kn/Ln  .  The  yd(m,  co )  is  spatially  re¬ 
sampled  by  performing  the  following  steps. 
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1.  Insert  (Kn-1)  zeros  between  each  pair  of  measured  spatial  samples  giving 

{nr 

yd(—’COn)f°r-m  =  ^±Kn’±1Kn’---- 

0  otherwise 

2.  Filter  z  (m,  a>n )  using  a  low-pass  filter  with  FIR  hn(m)  designed  to  approximate  the 
frequency  response.  Convolve  zip(m,  con  )  with  hn(m)  yields  the  low-pass  filtered  output 
ziP(m,  con  ). 

3.  Decimate  zip(m,  con  )  by  a  factor  of  Ln  to  obtain  the  focused  spatial  data  sequence 
yd  (m,  0)n )  =  zip  (Lnm,  (On ) 

After  forming  the  spatially  re-sampled  Mxl  data  vectors 

Yd ( con )  =  [y<; (0, con ),yd(l,(On),..., yd (M  - 1,  con )f 

For  n=l a  single  focused  covariance  matrix  at  frequency  O)0 is  estimated  by 

n= 1 

The  wide-band  source  location  can  be  computed  using  MUSIC  operated  at  frequency  a>0 . 

Simulation  Method 

A  16  element  line  array  with  sensors  spaced  a  half  wavelength  apart  at  (dma  =  n  was 
used.  The  array  focused  to  a>0  -  co^  -  tt  /  2  where  the  Rayleigh  resolution  limit 
corresponds  to  sin(  <9)=0.25.  For  each  trial,  32  independent  realization  of  yat  0)n )  were 
generated  for  each  O)n-Knco0/Ln.  Where  Kn=32  and  Ln=16,  17,...,  32.  In  these 
experiments  the  up-sampling  factor  Kn  factor  was  fixed  for  convenience,  which  led  to 
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con,  n=l,  17  non-uniformly  spaced  frequency  samples  spanning  the  band  from  n! 2 


to  K .  The  spectral  MUSIC  was  to  process  the  focused  covariance  matrix. 

Summary 

This  method  requires  interpolation  and  then  decimation  of  input  data  and  then 
computation  of  covariance  matrix.  The  interpolation  filter  was  designed  using  the 
minimum  mean  square  error  filter  design  method.  An  FIR  filter  of  length  2.J.  Kn  +1 
points  was  used  to  approximate  an  ideal  low  pass  filter  with  cutoff  frequency 
coc  -  j5n!  Kn .  P  is  between  0  and  1  (It  was  chosen  as  0.4).  The  filter  length  factor  J  is 

chosen  as  4.  The  computation  of  covariance  matrix  is  done  on  B  length  of  data  which  is 
done  across  various  frequencies.  So  the  value  of  B  could  be  large  in  the  case  of  wide¬ 
band  data.  The  MUSIC  algorithm  is  applied  to  resolve  various  DOAs.  Technique  looks 
good  some  work  need  to  be  done  to  resolve  value  of  B,  K  and  L  which  are  also  assumed. 
Selection  of  B,  K,  and  L  may  be  tricky  and  this  approach  then  will  not  be  useful  in  a 
generalized  case.  No  attempt  was  made  to  run  MATLAB  program  due  to  above  reasons. 
One  advantage  of  this  approach  is  that  it  does  not  require  preliminary  DOA  estimates  and 
iterations. 

2.5  New  Signal  Subspace  Direction  Of  Arrival  Estimator  for  Wide-band  Sources  by 
Yeo-Sun  Yoon,  Lance  M.  Kaplan  and  James  H.  McClellan  (Review)  [20] 

This  work  considers  a  uniform  linear  array  with  M  sensors  and  assumes  that  the 
frequency  bands  of  the  D  sources  are  known  [20,  38,  40],  They  are  also  overlapped.  The 
sensor  outputs  are  decomposed  either  using  a  filter  bank  or  FFT.  The  sensor  output  at 
frequency  rq  is  given  as 

x{a$  =  Ai(0)S{a$  +  ri{a$  (2-43) 
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The  relation  between  array  manifolds  of  different  frequencies  and  DO  As  is  as  follows: 


a(a>x,0x)  =  &(a>y,0y)a(a>z,0z)  (2.44) 

where  4>(ct)y,0y)  -  diag{e  jC0yTy,e  it0y2Ty i0,yMTy} 

The  relation  between  frequencies  and  DOAs  is  as  follows: 

=  coy  +  coz 

sin  0 '  =  —  sin  0  +  —  sin  6 


Let  X|=x(  cot )  where  cot  for  i=l,2,..K  is  within  the  frequency  bands  of  all  sources.  The 
number  of  frequency  bins  K  is  constrained  by 

K  >  maxi — — — ,3} 

M-D  J 


The  covariance  matrix  R,  can  be  defined  as: 

R,  =  E[xtxf  ]  =  AtRsi A?  +  a2 1  (2.45) 


where  Rs  i  =  E[s{coi)s{cof  ]  and/?  is  a  full  rank  matrix. 


There  are  D  largest  eigenvectors  corresponding  to  D  sources  and  their  range  is  the  same 
as  the  range  of  A,-.  Let  e  be  the  ordered  eigenvectors  of  /?,  from  the  largest  to  the  smallest. 
They  define  matrices  Fi  as  signal  range  space  and  W,  is  the  null  (noise)  space 


Fj 


—  [ei,D+Vei,D+2T"ei,M^ 


(2.46) 


(2.47) 


Then  following  can  be  specified  for  signal  and  null  range  spaces: 
Range  {F;}=  Range  {A,} 

Range  {lTi}=Null  Range  {A,H) 


43 


Let  A  (O  =  0)j  -  cOj  then, 


Range{<A>(Aco,  60)Fi }  =  Range{Aj(6)} 

where  0  =  [Ol,...QDf  and  dd  =  arcsin{— sin#rf  +  ^^sin#0} 

This  theorem  informs  that  a  signal  subspace  of  one  frequency  bin  can  be  linearly 

transformed  into  that  of  other  frequency  with  modified  DOAs  in 

,  „  Q)  .  n  a  .  n 
sin#  =^Lsin  #  H — -sin# 

If  #0  in  Range{<f>(Ao),  #0  )FI }  =  Range{Aj{6)}  is  the  same  as  one  of  the  0d  in  the  original 
signal  subspace,  this  DOA  is  preserved  in  the  new  signal  subspace. 

Assume  that  K  >  maxJ  — — — ,3 1  holds. 

[M-D  J 

Let  Ei  for  i=l,2. .  .K  be  Dx(M-D)  matrices  such  that 
E.  =  F^iAco^W, 


where  Aco  -cot-co^  Define  the  Dx(M-D)  matrix  B  such  that 

B=[E2  E3  ...  Ek] 

Then 


rank{B(6)B(0)H } 


D-\..JfO  =  Od 1 

D . if8*8d} 


Proof:  As  Range  {Wi}=Null  Range  { AiH}  has  been  specified  earlier  then 

adm)H%=  °r 


This  is  true  for  all  d.  We  also  know  that 
F^iA^Of  =TtHA(6>,)H 
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where  0.  d  -  arcsin<|  —  sin  0A  + 


u 


l 


CO. 


CO. \ 


sin# 


if  o  =  od  and  e2d=..=ekd=ed 
Therefore, 


Ei  =  Ti 


H 


aH  (a>i,dt)Wl 


aH  (a, 0d)Wi 


aH  {coi,6d)Wi 


(2.48) 


Giving  the  following  relation 

Et  =  T» 

This  is  the  dth  row  where  we  have  0T. 

Since  there  are  multiple  sources  then  there  is  a  possibility  that  one  of  the 

6  s  would  be  the  same  as  <9; .  If  we  use  atleast  three  frequency  bins  then  this  ambiguity 


* 

0T 

* 


can  be  removed.  Therefore  only  for  0  -  6d ,  B  becomes 


B  =  [E2..EK]  =  TlH 


and  it  loses  rank. 

The  estimation  process  is 


(2.49) 


45 


1.  Divide  the  sensor  output  into  J  identical  blocks 

2.  Compute  FFT  of  each  block 

3.  Find  xi  for  pre-selected 

4.  Find  the  signal  subspace  Fj  and  noise  subspace  Wi  by  eigendecomposition  of  the 
covariance  matrix  Ri 

5.  Find  E,  using  Ej  =  Ft"(p(Aa)i,0)HWi 

6.  Find  0  such  that 

0  =  arg  max  K{B(0)B(0)H } 

Summary 

This  method  requires  that  preliminary  estimate  of  the  DOA  that  could  be  one  of  the 
angles  in  the  estimation  and  number  of  sources  need  also  be  estimated.  There  is  missing 
information  in  this  work  so  no  further  study  is  planned  at  this  time  for  this  work.  We  may 
come  back  later  on  and  have  a  second  look  at  it  as  things  progress  with  other  work. 
Authors  have  extended  their  work  for  arbitrary  shaped  multidimensional  arrays.  They 
have  also  named  this  technique  as  Test  of  Orthogonality  of  Projected  Spaces  (TOPS)  [20, 
38-40], 

2.  6  A  Method  for  Wide-band  Direction  of  Arrival  Estimation  Using  Frequency- 
Domain  Frequency-Invariant  Beamformers  by  Tuan  Do-Hong,  Franz  Demmel, 
Peter  Russer  (Review)  [21-25] 

A  new  method  for  wide-band  DOA  estimation  using  arbitrary  antenna  array  based 
on  Frequency-Domain  Frequency-Invariant  Beamformers  (FDFIB)  is  proposed  by  these 
authors  [21-25].  Earlier  a  beam-space  processing  using  Time-Domain  Frequency- 
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Invariant  Beamformers  (TDFIB)  was  proposed  [22],  This  method  does  not  require 
preliminary  DOA  estimates  and  has  less  computational  complexity  than  beamspace  CSS 
[19].  However  the  frequency  invariant  characteristics  depend  on  the  number  of  antenna 
elements  within  the  arrays,  on  the  geometry  of  arrays  and  on  the  design  of  filters  in 
TDFIB. 

In  this  work,  Frequency-Domain  Beamformers  (FDBs)  are  used  with 
appropriately  designed  weights  at  different  frequencies.  This  ensures  that  beam-patterns 
of  FDBs  remain  constant  over  the  frequency  band.  This  approach  is  then  termed  as 
frequency-invariant  beamformer  as  beam  patterns  are  independent  of  the  frequency.  This 
technique  transforms  the  element-space  into  the  beam-space  and  acts  as  spatial  processor. 
DOA  is  then  estimated  using  the  well  known  narrow-band  MUSIC  that  is  applied  in 
beam-space  at  single  selected  frequency.  The  selected  frequency  should  also  be  within 
the  bandwidth. 

Frequency-Domain  Frequency-Invariant  Beamformers  (FDFIB) 

An  array  of  M  identical  elements  and  a  beamforming  network  of  J  frequency 
domain  beamformers  are  considered  and  it  is  assumed  that  there  are  D  wide-band  sources 
located  in  the  far  field.  The  jth  beam  pattern  is  given  by: 

Bj{(ok)  =  wIj{(Ok)b{cok,Qlb)  (2.50) 

where  b(G\,CLh)  is  the  steering  vector  and  Q.b  =  [cos </>b  sin 6b, sin </>b  sin 6b, cos 0b\  with 
azimuth  and  elevation  at  the  bth  direction.  Azimuth  is  between  -;rto  n  and  elevation  is 
between  0  and  nil.  w''((Ok  )  is  the  weighting  vector  with  weight  at  frequency  o\  of  the 
jth  beamformer  at  mth  antenna  element. 
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The  weights  at  each  frequency  are  chosen  such  that  Bj(o)k )  =  Bj(co0) .  Where  (O0 
is  the  focusing  frequency?  The  weights  for  frequency-invariant  beamformer  at  each 


frequency  are  determined  as 


(0  (co,)  =  a  e 

jm  V  k '  m 


-j^d‘mab+j^-d‘m£ib 


(2.51) 


Where  am  is  the  amplitude  weighting  coefficients  and  its  value  as  unity  is  suggested  in 
this  work? 

Wide-band  Signal  Model  in  Element  and  Beam-Space 

Denote  x  (n)  as  the  output  of  the  sensor  arrays  which  is  Mxl  vector.  Sensor  array  output 
in  frequency  domain  can  be  written  as: 

X(cok)  =  ,  r d  )Sd  (cok )  +  N(cok )  (2.52) 

d= 1 

where  Tp  -  [cos (f>d  sin  6d , sin  (f)d  sin  6d , cos 0d  ]'  and  in  matrix  notation  above  equation 
becomes: 

X(cok)  =  A(cok)S(cok)  +  N(o)k)  (2.53) 

where  A  is  MxD  source  direction  matrix,  S  is  the  Dxl  vector  of  signals  at  inputs  of  the 
array  and  N  is  the  noise  matrix. 

A  J-beamforming  network  is  used  and  it  is  assumed  that  D  is  less  than  equal  to  J 
and  J  is  less  than  equal  to  M.  Output  of  J  beamforming  network  in  frequency  domain  can 
be  written  as: 

Y(cok)  =  [Yl(cok),..Yj{(Ok),.YJ{cok)]T  =CH  (cok)X(cok)  =  Ac(cok)S(cok)  + Nc(o)k)  (2.54) 
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where  Y  is  Fourier  coefficient  of  beam-space  signal.  Ac(a)k )  =  C" (o\  )A(o\ )  is  source 
direction  matrix  in  beam  space,  C(a>k) -[wl(a>k),...wj(cok),..wJ(cok)]T is  called  the 
beamforming  matrix  and  w  is  weighting  vector  of  jth  frequency-invariant  beamformer. 


The  source  direction  matrix  in  beam-space  is  constant  for  all  frequencies  within 
the  signal  bandwidth  as  this  beamformer  is  designed  to  be  frequency-invariant.  A  single 
direction  matrix  Ac  ( a>0 )  can  characterize  the  DOA  for  the  wide-band  case.  It  is 

customary  to  assume  that  the  signals  and  the  noise  are  uncorrelated.  The  cross  spectral 
density  matrix  in  beam-space  is  given  by: 

Ry  (cok )  =  E{Y(o)k  )Yh  (o)k ) }  =  CH  (cok  )RX  (0)k  )C{cok ) 


where  Rx{cok)  -  E{X(cok)XH (cok)} is  the  cross-spectral  density  matrix  in  element-space. 
The  wide-band  covariance  matrix  in  beam-space  can  then  be  written  as: 

Ry  =X^r(fi\)  =  Z  CH {(Ok)Rx{(Ok)C{(Ok)  (2.55) 

k=l  k=l 


A  dense  grid  of  I  angle  points  (or  spatial  frequency  points)  of  azimuth  and  elevation  is 
defined.  For  the  MUSIC  algorithm  the  DOAs  are  determined  by  searching  the  peak 
positions  of  the  spatial  spectrum. 


9  FDFIB -MUSIC 


(T,)  =  - 


flc(I^)flc(I^.) 
ac  (f)  )U NU %  ac  (T) ) 


(2.56) 


Where  U  is  the  noise  subspace  matrix  and  is  obtained  from  the  eigenvalue  computation  of 
R  matrix  and  a'('  is  the  Jxl  source  direction  vector  in  beam  space. 


Summary  of  the  algorithm: 

Compute  the  FFT  of  the  data  to  go  to  frequency  domain 
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Frequency  domain  data  is  fed  to  J  beamformers 

Compute  output  of  the  J  beamformers  in  parallel  using  steering  vector  and  weights 
Compute  covariance  matrix  from  the  outputs  of  J  beamformers 

Compute  MUSIC  algorithm  using  a  selected  frequency  on  a  dense  grid  of  I  angle  points. 
The  coordinates  of  the  array  (dm)  are  not  known.  Azimuth  estimation  using  a  uniform 
circular  array  of  9  elements  is  used.  The  co-ordinates  of  the  elements  are  normalized  over 
A  —  d  fh  and  the  radius  of  array  r=  A .  Three  uncorrelated  wide-band  sources  with 

normalized  frequencies  are  considered.  The  focusing  frequency  6)0  is  selected  at  0)h . 

2.7  Wide-band  Direction  of  Arrival  Estimation  and  Beamforming  for  Smart 
Antennas  System  by  Tuan  Do-Hong,  Peter  Russer  (Review)  [21-25] 

This  paper  is  similar  to  the  previous  paper  from  these  authors  [21]  and  their  two 
other  papers  [24-25]  with  following  differences.  It  assumes  an  arbitrary  Uniform  Linear 
Array.  It  uses  a  wide-band  beamforming  method  with  prescribed  narrow  main-beam 
width  and  Low  Sidelobe  Level  (SLL)  using  spatial  interpolation.  It  uses  a  spatial 
interpolation  process  consisting  of  two  FDFIBs.  The  first  FDFIB  is  based  on  a 
(prototype)  FDFIB  when  the  inter-element  spacing  of  d  is  replaced  by  Nd.  In  this  case  N 
is  an  integer  referred  to  as  expansion  factor.  The  required  main-beam  width  can  be 
obtained  by  adjusting  the  number  N  in  the  combination  with  spatial  widows.  The  second 
FDFIB  is  applied  at  the  output  of  the  first  FDFIB  to  attenuate  grating  lobes,  which  appear 
due  to  changing  of  the  spacing.  The  attenuated  level  depends  on  the  required  SLL.  As  a 
result,  the  beam-pattern  with  narrow  main-beam  width  (due  to  larger  spacing)  and  low 
SLL  (due  to  the  attenuation  of  the  second  FDFIB)  are  simultaneously  obtained  without 
increasing  the  number  of  antenna  elements.  The  use  of  spatial  interpolation  provides  for 
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the  specification  of  main-beam  width  while  at  the  same  time  allowing  for  the 
specification  of  SLL,  one  can  reduce  number  of  antenna  elements  and  the  corresponding 
RF  modules,  A/D  converters  etc.,  while  still  retaining  same  main-beam  width  and  SLL  as 
traditional  beamforming  method  that  requires  larger  number  of  elements.  Moreover  due 
to  larger  inter-element  spacing,  the  mutual  coupling  between  antenna  elements  can  be 
eliminated. 

Summary 

Four  papers  were  reviewed  from  these  authors  [21-25]  and  they  are  similar  in  their  work 
with  notational  changes  in  equations.  That  makes  it  difficult  to  separate  their  work. 
Simulation  results  are  also  somewhat  similar.  They  have  not  provided  a  step  by  step 
mechanism  for  their  algorithm.  The  bright  part  of  their  work  is  that  it  does  not  require 
preliminary  estimate  of  the  DOA.  It  does  require  array  geometry  and  guessing  of 
common  frequency  so  their  algorithm  could  be  based  on  it. 

2.8  Theory  and  Design  of  Broadband  Sensor  Arrays  with  Frequency  Invariant 
Beam  Patterns  by  Darren  B.  Ward,  Rodney  A.  Kennedy,  Robert  C.  Williamson 
(Review)  [26] 

This  work  deals  with  the  problem  of  designing  a  uniformly  spaced  array  for  wide 
band  applications.  This  group  of  authors  has  published  three  papers  in  this  series  and  they 
will  be  reviewed  one  at  a  time  [26,29-30].  A  summary  of  all  these  work  will  be  provided 
at  the  end  of  the  third  review. 

Consider  broad  band  arrays  in  which  there  is  little  or  no  frequency  variation  in  the 
far-field  array  beam  pattern  over  an  arbitrarily  wide  desired  bandwidth.  The  asymptotic 
theory  of  unequally  spaced  arrays  is  used  to  derive  relationships  between  beam  pattern 
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properties  and  array  design.  Theory  in  this  work  assumes  planar  arrays  and  sources  are  in 
far  field.  Broadband  FI  array  is  defined  in  terms  of  the  array  beam  pattern.  These  beam 
patterns  are  assumed  to  be  identical  at  different  frequencies  and  require  a  compound 
array  of  k  subarrays.  Subarrays  are  identical  and  their  spatial  coordinates  are  expressed  in 
wavelength.  There  infinite  number  of  subarrays  would  be  required  for  producing  an 
identical  beam  for  wide  range  of  frequencies.  First  of  all  a  continuous  sensor  is  developed 
which  will  produce  frequency  invariant  beam  pattern  and  then  it  is  approximated  with 
group  of  discrete  sensors. 

First  consider  a  one-dimensional  (linear)  continuous  sensor  aligned  with  the  x-axis.  The 
output  of  this  continuous  sensor  is 

Zf  =  ]s(xJ)p{xJ)dx  (2.57) 

where  S  is  the  signal  received  and  p  is  the  sensitivity  distribution 

The  p  (x,f)  is  considered  as  the  aperture  distribution.  It  is  assumed  that  the 
sensitivity  distribution  is  absolutely  integrable  and  the  integral  exists  for  finite  power 
signals.  The  output  of  the  sensor  when  subject  to  plane  waves  arriving  from  an  angle  6 
is  given  by 

S(x,f)  =  e~i2m~'fxsiae 

where  c  is  the  speed  of  wave  propagation. 

The  output  of  the  sensor  is  a  function  of  6  and  the  sensor  beam  pattern  at 
frequency/  is  as  follows: 

bf(0)=  J  e-J2x~lfxsinSp(x,f)dx 
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A  wide-band  frequency  invariant  (FI)  sensor  will  have  pattern  as  frequently  invariant  and 
is  defined  as  ,  bf  (6)  =  b(Q) ,  for  all  />0. 

The  sensitivity  distribution  of  a  one-dimensional  sensor,  which  is  a  function  of  distance  x 
along  the  sensor  and  frequency/is  given  by 

P(x,f)  =  fG(xf  )  for  all  f>0 

where  G  is  an  arbitrary  absolutely  integrable  complex  function  of  a  single  real  variable. 
Then  the  far-field  beam  pattern  bf{0)  which  is  a  function  of  the  angle  #  measured 

relative  to  broadside  and  frequency/,  will  be  frequency  invariant. 

bf  (6)  =  b{0 )  =  |  e~i2m:^siDdG^)d^  (2.58) 

where  /  =xf. 

Above  mathematical  derivation  was  presented  as  a  theorem  by  authors.  A 
sensitivity  distribution  theorem  was  also  developed  by  authors.  Assume  b(0)  as  an 
arbitrary  continuous  square  integrable  frequency  invariant  far  field  beam  pattern  which  is 
specified  for  6e  (-nl2,nlT)  and  it  determines  p  (x,f)  uniquely.  Then  the  sensitivity 
distribution  p  (x,f)  of  a  linear  sensor  which  realizes  this  beam  pattern  must  satisfy  the 
following  conditions: 

1.  p(x,f)  -  fG(xf)  for  some  function  G. 

2.  G  has  a  Fourier  transform  F  satisfying 
T (s)=B(s)=b{sin  1  (sc)},  se  (-l/c,l/c) 

T(s)=A(s),  (-l/c,l/c) 

where  c  is  the  speed  of  wave  propagation  and  A{.)  is  an  arbitrary  square  integrable 
function  such  that 
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A[(-1)7c]=  lim  B(s ) 

.  ,(-!)' 
c 

for  i=0, 1 

Thus  the  only  freedom  in  choosing  p(x, f)  for  a  desired  FI  beam  pattern  is  in  the 
sufficiently  high  spatial  frequency  behavior  of  G.  b{6)  for  Oe  {-7112,7112)  determines 
p  (x,f)  uniquely. 

Assume  D-dimensional  continuous  sensor.  Let  the  output  of  the  sensor  be  given  by 

Z7  =jRS(x,f)p(x,f)dx 

The  sensor  has  a  frequency  invariant  far-field  beam  pattern  if 

p{x,f)  =  fDG{xf ) 

If  G  is  an  arbitrary  absolutely  integrable  complex-valued  function  then  it  can  be 
expressed  as: 

G(xf)=  Af{x)  =  Hx(f) 

A  defines  the  aperture  distribution  function  at  a  nominally  fixed  frequency/ and  H 
defines  the  primary  filter  at  a  single  point  x  on  the  sensor.  The  total  filtering  required  at  a 
fixed  point  x  can  be  expressed  as 

p{x,f)  =  fDHx(f ) 

where  fD  is  considered  as  a  secondary  filter.  It  is  independent  of  the  sensor  spatial 
vector  x  and  it  is  a  function  of  the  sensor  dimensions  D  only. 

G  is  a  symmetric  function  of  spatial  variable  x  and  of  the  frequency  variable  /. 
This  implies  that /and  x  can  be  interchanged  without  affecting  the  value  of  the  function. 
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Their  values  can  be  varied  one  at  a  time  and  keeping  the  other  constant.  The  primary 
filter  response  takes  the  same  shape  as  the  aperture  distribution. 

If  H  denotes  the  frequency  response  of  the  primary  filter  at  a  point  x  and  A  denotes  the 
aperture  distribution  for  a  given  frequency  then  H=A.  (fx)  (ignoring  some  mathematical 
notations). 

The  primary  filter  response  required  at  point  x  can  be  obtained  by  taking  a  slice 
through  the  aperture  distribution  from  the  origin  in  the  direction  of  x.  The  aperture 
distribution  can  be  determined  from  the  desired  beam  pattern  and  vice  versa.  The 
correspondence  between  aperture  distribution  and  primary  filter  response  is  for  both 
magnitude  and  phase.  All  primary  filter  responses  in  a  D-dimensional  frequency  invariant 
broadband  sensor  for  a  given  x  are  identical  up  to  frequency  dilation. 

If  we  concentrate  on  single  sided  one  dimensional  array  apertures  with  the  first 
element  located  at  x=0.  An  array  of  sensors  can  only  approximate  the  ideal  broadband 
continuous  sensor.  This  reduces  to  a  numerical  approximation  uniformly  in  /  to  the 
following  integral  representing  the  output  of  the  ideal  continuous  sensor  for  an  arbitrary 
signal  S. 

Zf  =  ]s(x,f)fG(xf)dx  (2.59) 

This  is  for/greater  than  0.  Assume  (x,  }  denote  a  finite  set  of  N  discrete  sensor  locations. 
In  approximating  a  finite  number  of  sensors  and  limiting  the  range  of  frequency  to  \fu  fu\ 
then  we  have: 

Zf=/'tg,SU,J)G{xJ)  (2.60) 

(=0 

S  is  the  complex  signal  received  at  point  x;. 
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G(x{f)  is  the  sampled  value  of  G(xf)  at  x  =x;- 

gi  is  a  frequency  independent  weighting  function  to  compensate  for  the  possibly  non- 
uniform  sensor  locations. 

An  important  aspect  of  this  broadband  array  design  is  that  the  array  design  comes  from 
approximating  an  integral  describing  a  broadband  FI  continuous  sensor.  A  trapezoid 
integration  method  is  used. 

Using  G(xf)=  Af(x)  =  Hx(f)  write  the  output  of  the  primary  filter  attached  to  the  ith 
sensor  as 

yi(f)  =  HxXf)S(xi,f) 

This  can  also  written  using  filter  dilation  theorem  as 

yAf)  =  nXi  (*,-/*!  f)S(.xt,f) 

This  emphasizes  that  only  one  primary  filter  shape  is  required  in  the  numerical 
integration  approximation  and  is  written  as 

Zf=fy(f)'Tx 

The  weighting  function  g  can  be  seen  to  relate  to  7x  via  an  un-illuminating  formula.  The 
weighting  functions  can  be  a  function  of  one  or  more  discrete  sensor  locations  but  are 
independent  of  the  frequency. 

It  is  assumed  that  the  aperture  distribution  is  a  slowly  varying  function  with  respect  to  x 

compared  to  the  exponential  term  in  bf(9 )  =  J  e~j2m:  fxsm0p(x,f)dx . 

The  block  diagram  of  general  single  sided  one  dimensional  broadband  FI  array  would 
look  like  a  FIR  block  diagram. 

1 .  The  primary  filters  are  simple  dilations  of  a  single  frequency  response. 


56 


2.  H(f)  =  HXi(f) 

3.  The  primary  filter  frequency  response  H(f)  is  similar  to  the  continuous  aperture 
distribution  shape  both  in  magnitude  and  phase. 

4.  The  primary  filter  outputs  can  be  combined  via  frequency  independent  weight  g 
that  depends  only  on  the  sensor  locations  generating  a  scalar  output. 

5.  All  sensors  share  a  common  secondary  filtering  response /to  generate  the  final 
output. 

2.9  FIR  Filter  Design  for  Frequency  Invariant  Beamformers  by  Darren  B.  Ward, 
Rodney  A.  Kennedy  and  Robert  C.  Williamson  (Review)  [29] 

This  work  uses  a  Frequency  Invariant  Beamformer  (FIB)  approach  and  in  this 
case  the  response  of  the  beamformer  is  constant  over  an  arbitrarily  wide  design 
bandwidth.  Previously  these  authors  have  presented  an  analog  technique  based  on 
approximating  an  ideal  continuous  aperture  [26].  This  filter  design  approach  just 
described  in  this  report.  Based  on  this  theory,  two  methods  namely  multirate  sampling 
and  single  sampling  rate  of  designing  FIR  filters  for  use  in  an  FIB  are  proposed  in  this 
paper  [26,  29-30], 

The  response  of  a  linear  continuous  aperture  to  planar  waves  from  an  angle  6 
measured  to  broadside  was  given  as: 


*^max 

mf)=  J 


jlirfx  sin  6 


p(x,f)dx ) 


where  p(x,  f )  is  the  aperture  illumination. 

It  is  a  continuous  function  of  both  location  x  and  frequency  /  and  c  is  the  speed  of 
wave  propagation.  The  response  remains  constant  if  the  aperture  illumination  is  given  by 
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f.G(xf),  where  G(.)  is  an  arbitrary  absolutely  integrable  function.  This  allows  breaking  of 
filtering  of  an  FIB  into  two  parts  namely  primary  filter  response  and  the  secondary  filter 
response. 

1 .  The  primary  filter  response  is  Hx(f)=G(xf) 

2.  The  secondary  filter  response/ is  independent  of  position. 

An  important  feature  of  the  FI  aperture  is  that  all  primary  filters  are  related  by 
dilation,  i.e.  if  Hx(f)  is  the  primary  filter  response  at  an  arbitrary  point  on  the  aperture, 
then  the  primary  filter  response  at  a  point  yx,y>0  is  given  by 

Hp(f)  =  G(ycf)  =  .  (2.61) 

The  continuous  aperture  is  approximated  using  numerical  approximation  of  the 
integral  described  in  the  above  equation  to  be  able  to  design  a  practical  frequency 
invariant  beamformer.  Let  xn  denote  a  set  of  N+l  sensor  locations  with  the  zeroth  sensor 
located  at  xo= 0.  If  we  limit  the  frequency  to  the  range  [//  fu],  we  can  use  the  following 
approximations: 

m  =  fi8nHn(f)ej2^  (2-62) 

n= 0 

where  Hd)  is  the  approximate  frequency  invariant  response 
g  is  a  spatial  weighting  term 

H  n  (/)  -  G(xnf)  is  the  primary  filter  response  of  the  nth  sensor. 

In  other  work  [26],  it  was  shown  how  to  obtain  g  for  the  case  corresponding  to  the 
trapezoidal  integration  method.  The  set  of  sensor  locations  can  be  determined  by 
minimizing  the  number  of  sensors  required  while  avoiding  spatial  aliasing. 

The  sensor  locations  are  given  by 
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x„  =  < 
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n, 


P  —  (  —  )n~P 

2  P-1 


(2.63) 


where  P  is  the  aperture  length  measured  in  half  wavelength.  Au  is  the  wavelength 


corresponding  to  the  upper  frequency  of  operation  and  N  =  P  + 


log 


(  P  ' 


Sl 


log 


\Jl  j 


vP-ly 


Design  of  the  primary  filters 

The  primary  filters  of  an  FIB  have  an  important  property  of  frequency  dilation.  It 
means  that  all  primary  filters  are  derived  from  a  single  reference  frequency  response,  and 
hence  all  primary  filter  coefficients  may  be  derived  from  a  single  set  of  coefficients. 
There  are  two  design  techniques  described  by  this  work  and  they  are  multirate  method 
and  single  rate  method. 

Multirate  Method 

A  primary  filter  response  Hre/f)  at  a  reference  location  xref  having  a  sampling 
period  T  will  have  href[k]  as  a  set  of  1  filter  coefficients.  The  primary  filters  needs  to  have 
the  required  dilation  property  if  the  nth  primary  filter  response  is  given  by 

(L- 1)12 

Hn(f)=  YuKeA^e-^  (2.64) 

k=-(L- 1)/2 

where  Tn=Txn/xref  is  the  sampling  period  of  the  nth  sensor. 

Multirate  sampling  is  generally  achieved  by  sampling  every  sensor  at  the  highest 
rate  required  and  then  use  decimation  to  bring  down  the  sampling  rate  to  the  desired 
sampling  rate.  Therefore  each  of  the  primary  filters  would  be  implemented  by 
•  Downsampling  by  yn  -  xn  /  xref 
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•  Applying  the  reference  primary  filter 

•  Upsampling  by  yn . 

The  aperture  length  is  defined  to  be  P  half-wavelengths  at  all  frequencies  within  the 
design  band.  Therefore,  the  nth  primary  filter  is  band  limited  with  Hn(f)=0  for 
|/|  >  Pc/( 2xn ) .  If  we  ignore  the  zeroth  primary  filter  (which  has  a  constant  response),  the 

primary  filter  with  the  widest  bandwidth  will  be  located  at  xl  -c  /(2fu ) .  It  will  have  the 

effective  bandwidth  as  Pfu,  requiring  a  sampling  rate  of  /i=2Pfu.  The  reference  primary 
filter  is  located  at  xKf=c/(2fu). 

Single  Rate  Method 

A  desired  primary  response  at  some  reference  location  will  have  a  set  of  reference 
coefficients  href[k]  .  The  nth  set  of  primary  filter  coefficients  can  be  obtained  by  first 
reconstructing  the  continuous  time  impulse  response  of  href[k],  applying  the  scaling 
property  of  the  Fourier  transform  and  resampling  the  scaled  impulse  response.  The 
primary  filter  coefficients  are  given  by: 

1  (i_1)/2  m 

hn[m]  =  —  X  KjWsmci - k) 

Yn  k=-(L- 1)/2  Yn 

where  yn  =  xn  /  xref  and  sin  c(v)  =  sin(^r)  Ktdc)  . 

The  reference  set  of  coefficients  must  first  be  convolved  with  the  coefficients  of  a 
low  pass  filter  having  a  cutoff  of  ynfs  1 2  to  avoid  temporal  aliasing.  For  yn  =  0,  the 
reference  coefficients  are  simply  an  impulse.  The  length  of  the  nth  primary  filter  should  be 
L  yn  a  predefined  number  of  coefficients  should  be  used  to  give  best  results. 
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If  the  input  signal  is  band-limited  to/jy,  the  minimum  sampling  rate  is  fs=2fu.  The 
location  of  the  reference  sensor  should  be  such  that  Hiej(f)=0.  Hence  for  the  single 
sampling  rate  FIB,  the  reference  sensor  is  located  at  xref=Pc/(2fu). 

Determining  the  Coefficients  of  the  Reference  Primary  Filter 

With  the  use  of  fG(xf)  and  the  change  of  variables  in  the  first  equation  in  this  work.  The 
Fourier  transform  relationship  between  the  desired  response  and  the  aperture  distribution 
is  apparent.  Thus  given  a  desired  response  and  the  visible  region  the  required  distribution 
is  given.  The  coefficient  of  the  primary  filter  can  be  calculated  with  a  vague  method 
described  by  the  author.  The  secondary  filter  is  a  differentiator  and  its  design  could  be  a 
Type  4  FIR  filter,  even  length  with  odd  symmetric  coefficients. 

Simulation  Example: 

Authors  used  an  aperture  length  of  P=4  half-wavelengths,  and  the  design  bandwidth  of 
200-3400  Hz,  requiring  17  sensors  and  a  total  array  size  of  3.4  m.  A  secondary  filter  with 
12  coefficients  and  a  uniform  aperture  illumination  were  used  in  both  cases.  They  show 
response  of  the  multirate  FIB  with  a  maximum  sampling  rate  of  30  kHz  and  a  reference 
filter  with  nine  coefficients.  The  single  sampling  rate  FIB  with  a  sampling  rate  of  8  kHz 
is  also  shown.  The  reference  filter  had  nine  coefficients  but  each  of  the  primary  filters 
had  a  minimum  of  51  coefficients  (with  the  larges  primary  filter  having  151  coefficients). 
Summary 

The  technique  presented  is  for  acoustic  wave  in  the  air.  This  approach  requires 
design  of  set  of  primary  filters  and  a  secondary  filter.  Design  process  is  described 
vaguely  and  its  implementation  in  MATLAB  would  be  a  difficult  task  as  of  missing 
information.  No  further  study  is  planned  for  this  work. 
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2.10  Broadband  DOA  Estimation  Using  Frequency  Invariant  Beamforming  by 
Darren  B.  Ward,  Zhi  Ding,  and  Rodney  A.  Kennedy  (Review)  [30] 

This  technique  is  for  wide-band  focusing  using  time  domain  processing.  It 
performs  beamspace  processing  using  frequency-invariant  beamformers  i.e.  beamformers 
whose  beampattems  are  constant  over  a  wide  frequency  band.  This  approach  uses  a  set  of 
appropriately  designed  beam-shaping  filters.  These  filters  ensure  that  the  same  array 
manifold  is  produced  for  all  frequencies  within  the  design  band.  The  proposed  estimator 
does  not  perform  frequency  decomposition.  This  approach  exploits  the  FIR  filtering 
based  approach  to  implicitly  perform  focusing  over  a  wide  frequency  band.  It  is  a  filter 
and  sum  approach.  Authors  have  provided  design  of  filters  for  this  scheme  in  their  earlier 
work  [26,  29]. 

Frequency  Invariant  Beamforming 

This  approach  uses  a  filter  and  sum  structure.  The  FIR  filter  at  the  mth  sensor  has 
response  as  Hm(f).  Hs(f)  is  an  optional  normalization  FIR  filter.  The  beam  shape  is 
constant  as  a  function  of  frequency  and  beam  shaping  is  performed  by  sensor  filters.  The 
response  of  this  beamformer  to  plane  waves  arriving  from  an  angle  6  is 

M 

r(0,  /)  =  //,(/)£  Hn  (2.65) 

m=l 

where  Tn(0)  is  the  propagation  delay  to  the  mth  sensor.  Above  equation  may  also  be 
written  as: 

r(Q,  f)  -  bH  (f)a(0,  /) 

Assume  that  each  row  of  b(f)  is  an  FIR  filter  with  filter  coefficients  bm[nc]  and 
nc=0,  .  NC-1.  These  filters  need  to  be  properly  design  so  the  response  of  the 
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beamformer  is  approximately  constant  with  respect  to  frequency  over  the  design 
bandwidth  giving  a  frequency  invariant  response. 

Frequency  Invariant  Beamspace  Processing  for  DOA  Estimation 
Consider  a  linear  array  of  M  sensors  that  are  not  necessarily  uniformly  spaced.  Assume 
that  D<M  far  field  broadband  signals  arriving  from  D  sources.  The  time  series  received  at 
the  mth  sensor  is 

y „  [*  1  =  Z \ Ik  - T.  (»)]  +  v. Ik]  (2.66) 

d-  1 

Its  frequency  response  can  be  written  as 

yj /]  =  fv  w-l"\[/]+v„[/]  (2.67) 

d= 1 

We  can  define  an  N-dimensional  vector  of  stacked  array  data  and  its  frequency  response 
can  be  written  as: 

y(f)  -  A(0,  f)s(f)  +  v(/)  (2.68) 

where  s(f)  is  the  Dxl  source  signal  vector.  A  is  the  MxD  source  direction  matrix  and  v  is 
the  noise  vector. 

The  source  signals  and  noise  have  finite  bandwidth  [fufu]-  We  want  to  determine 
the  source  direction  from  the  observed  array  data  vector  y[k]  over  a  finite  time  period 
n=l....N. 

Assume  we  apply  an  FIB  to  the  received  array  data.  The  beam  former  output  is 

M  NC- 1 

z[n]  =  X  X  bm\-nc- \ym\n-nc\  (2.69) 

m= 1  nc= 0 

where  b  is  the  set  of  FIR  filter  coefficients  on  the  mth  sensor.  The  frequency  response  of 
the  beamformer  output  is 
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z(f)  =  bH(f)y(f) 


(2.70) 


Assume  we  now  form  J  (D  <  J  <  M  )  such  beamformers  using  J  different  sets  of  filtering 
vectors.  Denote  the  stacked  vector  of  beamformer  response  as  z(f)  where  b  is  the  ith  set  of 
beam  shaping  filter  responses.  These  beamformers  are  designed  to  cover  a  spatial  vector 
in  which  the  sources  are  assumed  to  lie,  and  they  should  have  uniformly  low  side  lobes  to 
attenuate  unwanted  out-of-sector  sources  e.g.  Chebyshev  beamformers. 

Let  C(f)  be  the  MxJ  beamforming  matrix.  Because  the  beamformers  are  designed  to 
satisfy  the  frequency  invariant  property,  the  FIBs  source  direction  matrix  is 
approximately  constant  for  all  frequencies  with  the  design  band.  Hence  the  broadband 
source  directions  are  completely  characterized  by  a  single  beamspace  source  matrix  Ac. 
Assuming  the  source  signals  and  the  noise  are  uncorrelated  the  FIBS  data  covariance 
matrix  is 

Rz  (/)  =  E{z(f)zH  (/) }  =  A  (&)RS  ( f)Ac (0)"  +  Rv (/)  (2.71) 

The  broad  band  FIBS  data  covariance  matrix  is  now  in  a  form  in  which  conventional 
eigenbased  DOA  estimator  may  be  applied.  The  eigensubspaces  of  noise  and  signals  can 
be  obtained  and  MUSIC  algorithm  can  then  be  applied. 

Summary  of  the  proposed  algorithm 

Design  J  FIB’s  that  cover  the  selected  spatial  region. 

Calculate  the  broadband  FIBs  noise  covariance  matrix  Rv  (why  and  how  in  real  life  can 
this  be  calculated) 

Collect  data  from  each  of  the  J  beamformers  over  the  observation  period  n=l,...N  and 
estimate  the  broadband  FIBS  data  covariance  matrix  R?  using  the  following  relation. 
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Rz=^lLz[n]zH[n] 

N  n=  1 

where  z  is  computed  using  the  following  relation: 

M  NC- 1 

z[«]  =  XE  k\nch',M- 

m= 1  nc=0 

Using  Rz ,  we  can  form  signal  and  noise  subspace  and  use  MUSIC  algorithm  to  find 
DO  As. 

This  work  use  27  uniformly  spaced  elements  with  d  spacing.  Three  FIB’s  were 
designed  (using  FIR  filters  with  201  taps)  to  be  frequency  invariant  over  the  normalized 
frequency  band  [0.2,  0.4]  and  to  cover  the  spatial  sector  {80  to  100  degrees}.  The 
corresponding  beampattem  with  center  frequency  of  0.3  were  calculated.  Authors  show 
results  and  claim  that  they  are  better  than  CSS  and  Lee’s  paper. 

Summary 

This  approach  looks  very  attractive  and  needs  some  more  work.  It  is  still  not  clear 
as  how  to  design  these  FIR  filters  which  require  more  calculation  and  experimentation. 
We  need  to  find  a  way  to  design  these  FIR  filter  coefficients  need  to  be  more  focused  and 
find  a  way  to  compute  them. 

2.11  Cyclostationarity  Based  Coherent  Methods  for  Wide-band-Signal  Source 
Location  by  Giacinto  Gelli  and  Luciano  Izzo  (Review)  [39] 

This  work  exploits  cyclostationarity  properties  existing  in  modulated  signals  and 
uses  for  the  computation  of  direction  of  arrival  and  spatial  filtering  in  congested  areas. 
Signals  which  do  not  have  similar  cyclostationary  properties  do  not  pose  problems  for 
cyclic  methods  [39].  These  methods  can  also  be  used  where  number  of  sources  are 
greater  than  or  equal  to  the  number  of  sensors.  This  work  extends  the  cyclic  approach  to 
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wide-band  case  and  uses  Cyclic  Spectral  Density  Matrix  (CSDM)  for  the  estimating 
DO  As  at  different  frequency  values. 

Cyclic  wide-band  signal  source  location  method  should  have  following  properties: 

1 .  Exhibit  high  spatial  resolution 

2.  Suitable  to  work  in  a  multiple  access  scenario 

3.  Capable  of  also  working  in  a  fully  correlated  environment 

4.  Can  potentially  exploit  all  the  Signal  of  Interest  (SOI)  bandwidths. 

This  work  is  based  on  the  concept  of  focusing  transformations  [2,  20-21],  Different 
array  CSDMs  are  coherently  combined  into  a  single  matrix  using  a  common  frequency. 
This  is  done  by  means  of  appropriate  linear  transformations.  The  resulting  matrix  will 
condense  all  the  information  contained  in  the  frequency  dependent  array  CSDMs.  Then 
its  signal  subspace  properties  are  used  to  obtain  high  resolution  estimates  of  the  DO  As. 
The  scope  of  this  paper  is  two  fold: 

1.  To  show  theoretically  how  the  coherent  or  focusing  approach  can  be  applied  to 
the  cyclic  case. 

2.  To  extend  practical  focusing  techniques  from  the  conventional  to  the  cyclic  case. 
Two  techniques  are  proposed  in  this  paper. 

1.  It  is  an  extension  of  the  Wang  &  Kaveh’s  CSS  method  [6-7].  It  is  called  Cyclic 
Coherent  Signal  Subspace  (CCSS)  method. 

2.  The  second  part  uses  a  class  of  spatial  resampling  method  and  it  is  an  extension  of 
the  Array  Manifold  Interpolation  (AMI)  method  of  Krolik  [17-18].  It  is  renamed 
as  Cyclic  AMI  (CAMI) 


66 


Both  these  methods  can  perform  high-resolution  wide-band  signal  source  location 
even  in  fully  correlated  environment.  They  exploit  all  the  bandwidth  of  the  SOIs  and  are 
able  to  separate  SOIs  with  similar  characteristics. 

Background  on  Frequency  Domain  Cyclostationarity 

Second  order  wide  sense  joint  characterization  of  two  scalar  zero  mean  discrete 
time  complex  valued  signals  requires  knowledge  of  the  cross  correlations. 

Rxlx  2  (ft,  m)  =  E[xx  (ft  +  m)x*2  (ft)] 

Rxlx2  (ft,  m)  =  E[xx  (ft  +  m)x2  (ft)] 
where  E  denotes  statistical  averaging. 

If  xj  and  X2  are  jointly  second  order  cyclostationary  then  RxiX2  (n,m)  is  a  periodic  function 
of  n.  It  then  allows  a  generalized  Fourier  series  representation  with  respect  to  n  for  every 
value  of  m. 

a 

where  each  Fourier  coefficient 

R;uv(m)  =  {R,u„(n,m)e-*™) 

It  denotes  infinite  time  averaging,  represents  the  cyclic  cross  correlation  function 
and  the  values  of  a  are  referred  to  as  the  cyclic  frequencies.  When  xi(n)=X2(n),  the 
above  functions  R  reduce  to  the  Cyclic  Autocorrelation  Function  (CAF)  and  the  Cyclic 
Conjugate  Correlation  Function  (CCCF).  The  Fourier  transform  S  of  R  is  referred  to  as 
the  cyclic  cross  spectral  density  function  and  cyclic  conjugate  cross  spectral  density 
function.  When  a=  0,  both  R  and  S  will  reduce  to  the  conventional  cross  correlation 
function  and  the  power  cross  spectral  density  functions  respectively. 
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If  we  assume  that  xj(n)  and  X2< n )  satisfy  appropriate  mixing  conditions,  the 
functions  S  can  be  given  an  alternative  interpretation  in  terms  of  spectral  correlation. 
Then  Xi(v)n  and  X2(v)n  represents  their  Fourier  transform  of  finite  sequences.  It  results 
that 

<y)NxF  (H(v  -a))N]  (2.72) 

S  is  interpreted  as  the  limit  spectral  correlation  between  pairs  of  spectral 
components  of  xj(n)  and  x2(n)  spaced  a  cycles  apart.  Such  a  correlation  is  zero  for  a 
wide-sense  stationary  signal. 

Let  an  array  of  M  sensors  receiving  D  signals  and  D<M  which  are  bandlimited 
zero  mean  SOIs.  It  is  assumed  that  sources  are  in  far  field  and  wavefronts  are  planar. 
They  have  common  frequency  spectral  support  and  exhibit  cyclostationarity  with  a 
common  cycle  frequency.  Continuous  time  signal  received  at  the  mth  array  sensor  will 
have  the  following  form. 

*.(»)  =  IX  )]+<;,(')  (2-73) 

d= 1 

It  is  assumed  that  sensor  outputs  are  first  down  converted  to  baseband  and  successively 
sampled  and  processed  digitally.  Denoting /S=1/TS  >2W  the  sampling  frequency.  Then 
the  discrete  time  complex  envelope  at  the  mth  sensor  is  expressed  as: 

x,  (n)  =  fJsd[n~mk  (0d  )]^'2“(<W)  +  ik  in)  (2.74) 

d=\ 

where  v0  =  f0  /  f  s  and  mk  (6d  )  =  Tk  (dd  )/Ts  .  We  also  assume  that  m  term  is  integer. 

Moreover  the  discrete  time  SOIs  will  exhibit  cyclostationarity  with  discrete  time  cycle 
frequency 
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a  .  It  turns  out  that 


X(v)w  *  A(0,v  +  vo)S(v)„  +/(v)w  (2.75) 

which  are  Fourier  transform  and  expressed  in  frequency  domain.  If  we  assume  that  the 
interfering  and  noise  signals  are  uncorrelated  with  the  SOIs  and  do  not  exhibit 
cyclostationarity  with  the  considered  cycle  frequency  a ,  the  array  CSDM  defined  by  S  is 
given  by 

5 (v)  =  A(0,  V  +  v0  )S  (v)AH[*]  (0,  [-]  (v  -  or)  +  v0 )  (2.76) 

where  S  is  the  CSDM  of  the  SOIs  and  we  have  accounted  for  the  Hermitian  property  of  A 
with  respect  to  v.  Note  that  for  a  =0.  The  CSDM  given  by  above  equation  reduces  to  the 
conventional  array  spectral  density  matrix  (SDM).  Wide-band  cyclic  methods  which  are 
based  on  this  above  equation  are  expected  to  perform  well  even  in  unknown  or  time 
varying  interference  environments. 

If  we  assume  that  there  exists  at  least  a  value  of  v  such  that  the  array  manifold  is 
known  and  unambiguous  and  that  the  SOI  CSDM  has  full  rank  D,  the  signal  subspace 
approach  can  be  applied  to  obtain  high  resolution  estimates  of  the  signal  DOAs.  The 
wide-band  cyclic  MUSIC  method  estimates  the  signals  DOAs  by  resorting  to  the  SVD 
decomposition  of  the  array  CSDM  evaluated  at  a  single  value  of  v. 

Coherent  Cyclic  Methods  for  DOA  Estimation 

Authors  propose  to  perform  a  coherent  combination  of  the  array  CSDMs 
evaluated  in  correspondence  of  J  distinct  frequency  values  vj,  j=l,. . .  J.  Let 

Y(v)n=T(v)X(v)n 

This  represents  a  linear  transformation  of  X  and  T  is  a  nonsingular  MxM  matrix.  T 
should  satisfy  the  following  focusing  condition  [6-7]: 
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T  (v)A(0,  v  +  v0 )  =  A(0,  vf  +v0) 


(2.77) 


Where  v/  is  a  suitable  focusing  frequency  belonging  to  referred  focusing  bandwidth 
which  is  symmetric  to  v=0.  Denoting  y(n)  as  the  inverse  Fourier  transform  of  Y(v)  the 
array  CSDM  of  y(n)  is  given  by 

5“m(v)  =  lim  ^-£'[F(v)A,FH[*1([-](v-cr))/v] 

77L  J  N—>°°  Tv 

5^](v)  =  r(v)^T7£[x(v)wx"[*1(H(v-flf)^]*7,"[*,(H(v-ar)) 

77L  J  N-> oo  jy 

(v)  =  T(v)5“„  (v)  *  THn  ([-](v  -  a)) 

Modifying  this  by  using  Equation  12 

(v)  =  T(v)A(0,  v  +  v0 )S (v) Am*]  ( 0 , [-] {v  -  a)  +  v0 )T "[*]  (0,vf+vo)  (2.78) 

5*„(v)  =  A(0,vf  +  v0 )S"n  (v)Ah[*]  (0,  v f  +v0)  (2.79) 

Note  that  the  array  CSDM  at  frequency  v  has  been  focused  to  frequency  v/  by  means  of 
the  transformation  T(v).  Therefore  all  the  different  frequency  contributions  can  be 
coherently  combined  to  obtain  the  matrix 

Kn  =  (2.80) 

where  g(v)  is  a  complex  weight  function  to  be  suitably  chosen  and  the  discrete  sum 
ranges  over  J  frequency  values  belonging  to  Cl  f  .  By  substituting  (2.79)  into  (2.80)  we 

get  the  following: 

Km  =  M0,v,  +v„)R‘„MAHl‘\e,V/  +v0) 

where  R  is  the  DxD  matrix  that  combines  all  the  focused  SOI  contributions. 
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If  the  array  manifold  is  known  and  unambiguous  for  each  v,  the  matrix  A  has  full 


rank  D.  Then  under  the  assumption  that  R  also  has  full  rank  D  which  must  be  assured  by 
an  appropriate  choice  of  g(v),  it  results  that  R  has  rank  D<M  and  hence  any  signal 
subspace  method  can  be  applied  to  obtain  high-resolution  estimates  of  the  DOAs. 


]&n=USV“  = 


This  denotes  the  SVD  decomposition  of  R.  where  U  and  V  are  unitary  matrices  and  5  is  a 
real  non-negative  diagonal  matrix.  Since  R  has  rank  D  and  accounting  for  (15)  it  turns  out 
that  AH  {6, v ,  +  v0)U0  =  0  and  hence  the  maxima  of  the  spatial  spectrum 
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It  can  be  utilized  to  determine  the  unknown  DOAs. 

Conventional  focusing  methods  are  able  to  work  in  a  multipath  environment.  This 
property  is  also  extended  to  the  cyclic  case.  Sss  is  singular  at  every  frequency,  the 
frequency  averaging  performed  in  (16)  to  obtain  Rss  removes  the  singularity.  If  we  denote 
Av  as  the  spacing  between  two  consecutive  frequency  values,  for  Av « 1  and 
Qf  ~  Owe  can  write: 


PSS[*]  —  a  Z^(v)5ss„(v)Av 

AV  veil, 


‘ff 


P-ssi*]  ~  ^  |  <?^v)5ss[.:]  (v)Av 


Hence  choosing  g(v)  =  e  j2m'"'" ,  we  get 
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which  has  rank  D  and  mo  is  suitably  chosen  according  to  the  signal  characteristics.  Since 
this  approach  is  a  sub  space  based  approach  and  it  requires  estimation  of  D.  Generally 
this  can  be  done  using  AIC  or  MDL  algorithms  but  there  is  an  absence  of  noise  term  in 
Rxx  and  therefore  other  techniques  will  be  required. 

Practical  Focusing  Strategies 

Authors  suggest  that  T  should  be  estimated  either  using  the  Wang  &  Kaveh’s 
approach  [6]  for  single  group  case  or  use  Hung  &  Kaveh’s  approach  [7-8]for  multigroup 
case.  In  both  cases  preliminary  estimates  of  DOA  are  required. 

A  different  approach  which  does  not  require  preliminary  DOA  estimate  is  the 
class  of  spatial  resampling  method  of  Krolik  and  Swingler  [17-18].  Their  technique 
exploits  the  separability  of  the  array  manifold  with  respect  to  the  unknown  DOA  and  the 
frequency.  By  exploiting  the  series  of  expansion  of  a  plane  wave  in  polar  coordinates,  it 
can  be  shown  [32]  that  the  kth  element  of  the  steering  vector  {a(6d,v  +  v0)} in 

X(v)N  ~  A(0,v  +  v0)S(v)N  +  I(v)N  can  be  expressed  as: 

{a(0d,v  +  vo)}k  =  e^^)fM'c)  cos(4-*„) 

{a(0d,v  +  vo)}k  =  i'"7„[2j(v  +  v0)/J-]xe_Me;^ 

c 

where  (rk,(j)k )  are  polar  coordinates  of  the  kth  sensor  and  J„  is  the  nth  order  Bessel  function 

of  the  first  kind.  The  function  J  decays  faster  than  exponentially,  we  can  truncate  the 
infinite  sum  as: 

{a(ed,v  +  vQ)}k~  X  jnJn [2;r(v  +  v0 )fs  — ] x e~jn^ejn0d  (2.81) 

c 

n=-n£  ^ 
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which  allows  one  to  approximately  express  the  steering  vector  {a(0d ,  v  +  v0)}  in  separable 
form  as: 

{ a(0d,v  +  vo)}~G(v  +  vo)w(0d )  (2.82) 

where  G  is  an  Mx2Ne  matrix  whose  (k,n)th  element  is 

[G(v  +  v0)}kn  =  jnJn[2x(v  +  v0)fs  for  k  =1,..M, 

c 

where  as  w  is  the  2ne  column  vector  whose  nth  element  is  {w(0d)}n  —  e'"6’-' 

Therefore  by  substitution  of  (2.82),  the  focusing  condition  reduces  to 

T  (v)G(v  +  v0 )  =  G(v  f  +v0) 

Now  T(v)  can  be  solved  without  requiring  any  preliminary  DOA  estimate  using  G 
matrices.  This  method  is  referred  as  CAMI. 

Note:  Number  of  terms  in  summation  increases  with  increasing  size  of  the  array  and 
increasing  focusing  bandwidth. 

1.  For  a  fixed  number  of  M  it  puts  a  limitation  on  the  maximum  number  of 
terms. 

2.  The  accuracy  will  degrade  as  the  size  of  the  array  is  increased. 

Since  focusing  transformations  basically  exploit  the  spatial  properties  of  the  array,  it  is 
clear  that  almost  every  focusing  technique  can  be  extended  with  minor  modification  to 
the  cyclic  one. 

Summary 

This  approach  is  based  on  previously  reviewed  two  approaches  of  Wang  &  Kaveh  [2,20- 
21]  and  Krolik  [17].  They  proposed  ways  to  compute  DOA  using  above  mentioned 
approach.  Authors  specify  another  focusing  matrix  for  use  in  their  algorithm.  One 
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drawback  of  their  work  is  that  they  are  restricting  their  signal  to  cyclostationary  signals 
which  is  not  a  good  assumption.  We  are  looking  more  a  general  approach  and  hence  we 
are  not  pursuing  this  approach  at  this  time. 

2.12  Fabrizio  Sellone, ’’Robust  Auto-Focusing  Wide-band  DOA  Estimation  (Review) 
[31] 

Authors  propose  a  new  way  of  designing  focusing  matrices  which  has  some  robustness.  It 
does  not  need  initial  estimates  of  DOA.  He  also  claims  that  computational  requirement  is 
also  reduced  when  compared  to  RSS  and  SST  approach  [32,33]. 
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Chapter  3 


Simulation 

Simulation  is  the  best  way  to  verify  the  validity  of  any  algorithm.  Simulation 
program  first  of  all  requires  generation  of  appropriate  data  which  is  reasonable  for  the 
algorithm  under  the  test.  The  input  data  should  also  have  any  assumption  or  initial 
condition  that  are  necessary  and  should  be  spelled  out.  Algorithm  then  needs  to  be  coded 
properly  and  debugged.  Obtained  results  should  be  analyzed  to  check  validity.  Results 
can  then  be  shown  in  tabular  form  and  or  graphically  in  one  dimensional  or  three 
dimensional  forms.  Plotting  of  results  also  requires  appropriate  plotting  mechanism.  In 
this  work  it  is  best  to  test  our  algorithm  in  MATLAB  which  provides  mechanism  to 
generate  data,  commands  to  execute  any  algorithm  at  various  levels.  It  also  has  various 
ways  to  plot  results  from  one  dimensional  to  three  dimensional  plots. 


3.1  DO  A  estimation  for  narrow-band  sources 

A  Uniform  Linear  Array  of  sixteen  equally  spaced  Omni-directional  sensors  was 
used  as  shown  in  Figure  3.1. 


Figure  3.1  :  Uniform  Linear  Array  of  16  elements 
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The  spacing  between  sensors  is - 

^ft 0 


where  c  is  the  velocity  of  propagation  and 


f0  is  the  central  frequency.  The  signal  sample  size  is  4096.  The  MATLAB  simulation 

program  first  of  all  generates  random  data  using  two  sources  located  at  9  and  12  degrees. 
It  assumes  16  sensors  located  at  equidistance  in  an  Uniform  Linear  Array.  Gaussian  noise 
has  been  added  to  the  signals.  The  signal  to  noise  ratio  is  10  dB.  MATLAB  simulation 
using  above  mentioned  data  generation  was  performed  Figure  3.2  shows  DO  A  for  two 
sources  in  narrow-band. 
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Figure  3.2:  DOA  estimation  for  two  narrow-band  sources 
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3.2  DOA  estimation  for  wide-band  sources 


There  are  number  of  wide-band  DOA  algorithms  that  have  been  discussed  in  the 
previous  chapter.  Most  popular  wide-band  DOA  algorithms  are  signal  subspace  based 
algorithms.  These  subspace  algorithms  can  further  be  divided  into  incoherent  and 
coherent  signal.  MATLAB  programs  have  been  written  for  both  these  approaches. 

Wang  and  Kaveh  [6,  7]  proposed  Coherent  Signal-Subspace  (CSS)  method  for 
detection  of  DOA  for  wide-band  sources.  This  technique  separates  the  wide  frequency 
band  into  narrow-band  components.  The  data  set  for  this  algorithm  is  divided  into  64 
segments  and  each  segment  contains  64  samples.  It  also  assumes  an  Uniform  Linear 
Array  of  16  sensors.  The  time  domain  samples  will  be  transformed  into  frequency 
domain  by  applying  a  64  point  FFT  to  each  of  the  64  segments.  The  Coherent  Signal 
Subspace  approach  proposed  by  Wang  &  Kaveh  [6]  follows  following  computational 
steps. 

1 .  Compute  64  sets  of  64-point  FFT 

2.  Compute  33  Covariance  matrices  (16  by  16) 

3.  Computation  of  initial  DOA  estimate  using  MUSIC  algorithm  [6] 

4.  Computation  of  33  Focusing  matrices 

5.  Computation  of  Focus  matrix 

6.  Computation  of  number  of  sources 

7.  Separation  of  Signal  &  Noise  subspaces 

8.  Compute  DOA  using  MUSIC  algorithm 


In  order  to  demonstrate  the  performance  of  the  DOA  algorithm  for  wide-band 
signals,  an  Uniform  Linear  Array  of  sixteen  equally  spaced  Omni-directional  sensors  was 


used.  The  spacing  between  sensors  is  - ,  where  c  is  the  velocity  of  propagation  and 

Vo 
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fa  is  the  central  frequency.  Two  wide-band  sources  at  DOA  Ox  and  62  were  assumed. 

The  signals  are  stationary  zero  mean  band  pass  white  Gaussian  processes  with  central 
frequency  fa  =100 Hz  and  bandwidth  s  =  40m .  The  array  noise  is  also  stationary  zero  mean 
band  pass  with  the  same  pass  band  as  the  signal  with  a  SNR  of  KMB  at  each  sensor. 
Source  signals  and  the  noise  are  random  processes  with  a  bandwidth  of  40  Hz.  The 
sampling  frequency  is  chosen  to  be  300  Hz.  The  signal  will  be  observed  over  a  period  of 
T0  seconds  and  T  will  be  divided  into  k  =  64  segments.  On  each  of  those  segments,  the 

array  output  along  with  corresponding  noise  will  be  decomposed  into  narrow-band 
components  using  a  fast  Fourier  transform.  The  total  number  of  samples  taken  by  each 
sensor  will  be  4096. 

The  MATLAB  simulation  program  first  of  all  generates  random  data  using  two 
sources  located  at  9  and  12  degrees.  It  assumes  16  sensors  located  at  equidistance  in  an 
Uniform  Linear  Array.  Gaussian  noise  has  been  added  to  the  signals.  The  signal  to  noise 
ratio  is  10  dB.  The  data  set  consists  of  64  segments  of  64  sets  of  data  [6].  The  data  is 
bandpass  filtered  using  a  Butterworth  filter  and  has  wide-band  characteristics  of 
containing  frequencies  from  80  Hz  to  120  Hz.  MATLAB  simulation  using  above 
mentioned  data  generation  was  performed  Figure  3.3  shows  initial  estimate  of  DOA.  It 
can  be  seen  that  it  is  sort  of  pointing  towards  an  angle  of  10.5  degrees.  Using  this  initial 
estimate,  simulation  continues  to  perform  other  operations  and  was  able  to  show  two 
peaks  at  angles  of  9  and  12  degrees  very  clearly.  These  two  DOA  peaks  are  shown  in 
Figure  3.4. 
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Initial  Estimates  of  Wideband  DOA  MUSIC  Algorithm 


Figure  3.3:  Initial  estimate  of  wide-band  DOA  CSS  algorithm 
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Figure  3.4:  Final  estimate  of  wide-band  DOA  CSS  algorithm 
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Simulation  program  which  was  generated  earlier  for  the  Coherent  Signal 


Subspace  method  has  been  modified  to  incorporate  incoherent  signal  subspace  method 
proposed  by  Wax  [5].  The  simulation  program  uses  similar  input  data  as  was  used  in  the 
case  of  CSS.  Figure  3.5  shows  DOA  plots  showing  two  peaks  at  9  and  12  degrees  as  was 
obtained  earlier.  The  only  difficult  part  is  that  this  approach  uses  lot  more  computation  as 
MUSIC  algorithm  was  used  for  each  frequency  and  then  it  was  averaged. 
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Figure  3.5  :  DOA  estimate  using  incoherent  signal  subspace  method 

Simulation  program  which  was  generated  earlier  for  the  Coherent  Signal 
Subspace  method  has  been  modified  to  incorporate  bilinear  transformation  algorithm  as 
proposed  by  Shaw  [8  ].  The  simulation  program  uses  similar  input  data  as  was  used  in  the 
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case  of  CSS.  Figure  3.6  shows  DO  A  plots  showing  two  peaks  at  9  and  12  degrees  as  was 
obtained  earlier. 
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Figure  3.6:  DOA  estimation  for  wide-band  sources  using  Bilinear  Transformation 

Simulation  program  which  was  generated  earlier  for  the  Coherent  Signal 
Subspace  method  has  been  modified  to  incorporate  BASS-ALE  algorithm  as  proposed  by 
Buckley  [41].  The  simulation  program  uses  similar  input  data  as  was  used  in  the  case  of 
CSS.  Figure  3.7  and  3.8  shows  DOA  plots  showing  two  peaks  at  9  and  12  degrees  as  was 
obtained  earlier.  The  only  difficult  part  is  that  this  approach  uses  lot  more  computation  as 
MUSIC  algorithm  was  used  for  each  frequency  and  then  it  was  averaged. 
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Figure  3.7  :  DOA  estimation  for  wide-band  sources  using  BASS-ALE  algorithm 
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Wideband  DOA  using  BASS-ALE  Algorithm 
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Figure  3.8:  DOA  estimation  for  wide-band  sources  using  BASS-ALE  algorithm 

Simulation  program  which  was  generated  earlier  for  the  Coherent  Signal 
Subspace  method  has  been  modified  to  incorporate  wide-band  DOA  in  time  domain.  We 
inserted  16  delays  in  each  sensor  data  and  formed  a  data  vector  of  256.  This  then 
translated  into  a  covariance  matrix  of  size  256*256.  Eigendecomposition  was  again  on  a 
matrix  size  of  256  giving  256  eigenvalues.  MUSIC  algorithm  was  then  applied  to  get 
similar  DOA  plots  and  two  peaks. 
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Chapter  4 


Analysis  of  Computational  Requirements  of  Wide-band  DOA 

Algorithms 

Passive  detection  of  objects  has  been  studied  for  more  than  20  years  [1-14].  This 
approach  evades  detection  by  others  and  has  many  applications.  There  are  more  than 
thirty  publications  for  wide-band  detection  of  Direction  of  Arrival  (DOA)  algorithms 
which  are  available  in  the  literature.  These  algorithms  are  generally  presented  in  a  very 
complex  or  condense  form,  which  are  not  easily  understandable  for  people  who  are 
outside  that  narrow  field.  It  is  not  known  which  class  of  techniques  would  be  appropriate 
for  implementing  them  in  hardware  and  would  be  useful  for  real  time  applications.  One 
has  to  cut  through  all  the  mathematics  and  convert  algorithms  into  simple  arithmetic 
operations  before  architecture  can  be  visualized.  There  is  a  need  to  bridge  a  gap  between 
the  design  of  computer  hardware  especially  special  purpose  parallel  architectures  and 
available  algorithms  for  various  interdisciplinary  problems. 

This  work  is  the  first  step  in  sorting  out  which  algorithm  is  appropriate  for  further 
study  and  for  its  hardware  implementation  for  real  time  applications.  Wide-band  DOA 
algorithms  available  in  the  literature  have  been  reviewed  in  Chapters  1  and  2.  They  are 
also  simulated  in  MATLAB  and  their  results  are  provided  in  Chapter  3.  These  wide-band 
algorithms  can  be  divided  into  following  categories: 

•  Modal  Decomposition  Signal  Subspace  (In-coherent  signal  subspace) 

•  Coherent  Signal  Subspace  Method 

•  Rotational  Signal  Subspace  Method 
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•  Wide-band  DOA  detection  using  focusing  matrices 

•  Wide-band  DOA  detection  using  beamforming  approaches 

•  Combination  of  beamforming  and  focusing  approaches  restricted  to  chirp  or 
cyclostationary  signals. 

•  Use  of  ARMA  model  and  Bayesian  approaches 

•  Use  of  maximum  likelihood  algorithms 

•  Bilinear  Transformation  Method 

•  Wide-band  DOA  in  time  domain 


Signal  subspace  approaches  are  very  popular  for  computing  DOA  for  both 
narrow-band  and  wide-band  sources.  One  problem  with  them  is  that  they  may  not  be 
optimal  but  produce  computationally  efficient  algorithms.  Signal  subspace  based 
approaches  are  further  subdivided  into  coherent  based  and  incoherent  based  signal 
subspace  approaches.  Incoherent  signal  subspace  approach  decomposes  signals  into 
individual  narrow-band  frequencies  and  then  combined  them  to  produce  final  results. 
They  are  computationally  expensive  and  are  unable  to  resolve  coherent  sources  [2,  18- 
19].  Some  of  the  subspace  approaches  require  initial  estimates.  If  the  initial  estimate  is 
not  accurate  then  final  DOA  will  also  not  be  accurate  or  may  have  some  issues  with  it. 
There  will  also  be  issues  with  bias  and  variances. 

In  order  to  implement  DOA  algorithm  in  hardware  for  real  time  applications,  it  is 
important  to  use  a  computationally  efficient  algorithm.  One  approach  is  to  evaluate 
computational  requirements  of  currently  available  wide-band  DOA  algorithms  and  select 
one  of  them  for  hardware  implementation.  An  interdisciplinary  research  is  performed  for 
the  development  of  Parallel  Architectures  suitable  for  real  time  wide-band  digital  receiver 
applications. 
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Development  of  wide-band  DOA  algorithm  started  with  the  incoherent  signal 
subspace  approach  which  is  a  brute  force  extension  of  narrow-band  case.  This  approach 
works  but  computationally  expensive.  Wang  &  Kaveh  [6]  proposed  a  Coherent  Signal 
Subspace  approach  which  creates  a  focusing  matrix  using  a  single  frequency.  This  has  a 
simple  transformation  scheme  and  requires  initial  estimates  of  DOA.  Results  are 
reasonable  and  produce  two  peaks.  Hung  &  Kaveh  [7]  extended  their  previous  work  with 
a  promise  of  statistically  better  results.  They  introduce  concept  of  rotational  signal 
subspace  making  it  more  complex  or  accurate  focusing  matrix.  The  cost  to  their  approach 
is  that  they  added  a  computational  step  of  singular  value  decomposition.  There  may  be 
little  improvement  in  accuracy  of  results  with  additional  computational  cost.  Their  work 
was  simulated  in  MATLAB  and  results  were  obtained  and  were  shown  in  Chapter  3. 


Shaw  [8]  extended  work  of  Wang  &  Kaveh  [6]  and  Hung  &  Kaveh  [7]  and 
introduced  a  bilinear  transformation  approach  with  certain  approximation.  Advantage  of 
their  work  is  that  it  does  not  require  initial  estimate  of  DOA.  Their  work  was  simulated  in 
MATLAB  and  results  were  obtained.  Their  algorithm  is  very  sensitive  to  certain 
assumptions  and  parameter  values  which  make  it  unattractive  for  hardware 
implementation. 

The  algorithm  proposed  by  Ta-Sung  Lee  [19]  and  he  decomposes  data  using 
bandpass  filters  into  J  frequency  beams.  It  performs  beamspace  transformation  and 
computes  weights  using  least  square  method.  It  then  computes  beamspace  data  matrix 
and  focuses  on  single  reference  frequency  which  would  be  something  similar  to  CSSM 
method.  It  performs  transformation  into  K  beamspaces  and  forms  beamspace  data  matrix. 
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This  beamspace  data  matrix  then  focused  on  a  single  reference  frequency  out  of  J 
frequency  bands.  The  design  of  beamspace  data  matrix  is  described  which  requires  first 
design  of  beamforming  matrices.  They  are  again  designed  in  a  least  square  sense.  The 
problem  is  then  reduced  to  beamspace  data  correlation  and  noise  matrices.  Authors  then 
apply  their  own  derived  root  MUSIC  algorithm  which  could  be  substituted  with  the 
MUSIC  algorithm.  The  algorithm  does  not  require  any  preliminary  DOA  estimates. 

This  DOA  estimator  is  suboptimum  according  to  the  author.  It  would  also  result 
in  degradation  in  estimation  accuracy  at  low  SNR.  Even  with  the  true  DOA,  the  DOA 
estimates  may  not  be  exactly  identical  as  signal  roots  may  not  lay  on  the  unit  circle.  The 
proposed  parallelized  estimator  behaves  as  a  mixture  of  the  root-form  and  spectral-form 
estimator.  This  approach  provides  an  excellent  alternative  to  CSSM  but  frequency 
decomposition  and  calculations  of  weight  and  beamspace  data  matrices.  There  may  be 
some  alternative  ways  to  calculate  weight  and  compute  data  beamspace  matrices.  In  the 
end  it  again  uses  something  similar  to  MUSIC  algorithm  to  compute  DOA.  This 
algorithm  was  not  implemented  due  to  lack  of  information  and  we  are  looking  into  filling 
those  gaps.  The  previously  MATLAB  program  can  easily  be  extended  to  accommodate 
this  algorithm. 

We  have  investigated  work  of  Krolik  and  its  extension  [17-18].  Our  program  did 
not  give  two  peaks  as  expected.  Problems  could  be  in  missing  information  regarding 
transformation  matrix  that  could  lead  to  not  getting  appropriate  result.  Moreover  their 
paper  uses  random  processes  as  data  which  may  contribute  to  getting  two  peaks  in  their 
case  and  not  in  our  case.  This  technique  requires  computation  for  each  steering  direction 
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0  increasing  the  computation  requirements  but  eliminates  the  selection  of  initial  focusing 
angle. 

The  second  of  Krolik  work  uses  interpolation  and  then  decimation  of  input  data  and 
then  computation  of  covariance  matrix.  Technique  looks  good  some  work  need  to  be 
done  to  resolve  value  of  B,  K  and  L  which  are  also  assumed.  Selection  of  B,  K,  and  L 
may  be  tricky  and  this  approach  then  will  not  be  useful  in  a  generalized  case.  No  attempt 
was  made  to  run  MATLAB  program  due  to  above  reasons.  One  advantage  of  this 
approach  is  that  it  does  not  require  preliminary  DOA  estimates  and  iterations. 

The  work  of  Yoon  [20,  38-40]  requires  that  preliminary  estimate  of  the  DOA  that 
could  be  one  of  the  angles  in  the  estimation  and  number  of  sources  need  also  be 
estimated.  There  is  missing  information  in  this  work  so  no  further  study  is  planned  at  this 
time  for  this  work.  They  have  also  named  this  technique  as  Test  of  Orthogonality  of 
Projected  Spaces  (TOPS). 

Four  papers  were  reviewed  from  Tuan  Do-Hong  and  Peter  Russer  [21-25]  and 
they  are  similar  in  their  work  with  notational  changes  in  equations.  That  makes  it  difficult 
to  separate  their  work.  Simulation  results  are  also  somewhat  similar.  They  have  not 
provided  a  step  by  step  mechanism  for  their  algorithm.  The  bright  part  of  their  work  is 
that  it  does  not  require  preliminary  estimate  of  the  DOA.  It  does  require  array  geometry 
and  guessing  of  common  frequency  so  their  algorithm  could  be  based  on  it. 

Darren  Ward  et  al  [26,  29-30]  produced  a  group  of  three  papers  and  performs 
filter  and  sum  beamforming  in  frequency  invariant  fashion.  They  proposed  a  design  of 
FIR  filter  and  then  computed  covariance  matrices.  Three  FIB’s  were  designed  (using  FIR 
filters  with  201  taps)  to  be  frequency  invariant  over  the  normalized  frequency  band  [0.2, 
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0.4]  and  to  cover  the  spatial  sector  {80  to  100  degrees}.  The  corresponding  beampattern 
with  center  frequency  of  0.3  were  calculated.  This  approach  looks  very  attractive  and 
needs  some  more  work.  It  is  still  not  clear  as  how  to  design  these  FIR  filters  which 
require  more  calculation  and  experimentation.  Filter  specifications  were  not  provided. 
We  need  to  find  a  way  to  design  these  FIR  filter  coefficients  which  needs  to  be  more 
focused  and  find  a  way  to  compute  them. 

Gelli  and  Izzo  [39]  work  is  based  on  previously  reviewed  two  approaches  of 
Wang  &  Kaveh  [2,20-21]  and  Krolik  [17].  They  proposed  ways  to  compute  DO  A  using 
above  mentioned  approach.  Authors  specify  another  focusing  matrix  for  use  in  their 
algorithm.  One  drawback  of  their  work  is  that  they  are  restricting  their  signal  to 
cyclostationary  signals  which  is  not  a  good  assumption.  We  are  looking  more  a  general 
approach  and  hence  we  are  not  pursuing  this  approach  at  this  time. 

Sellone  [31]  proposed  a  new  way  of  designing  focusing  matrices  which  has  some 
robustness.  It  does  not  need  initial  estimates  of  DOA.  He  also  claims  that  computational 
requirement  is  also  reduced  when  compared  to  RSS  and  SST  approach  [32,33]. 

Wide-band  DOA  algorithms  reviewed  in  Chapter  2  and  simulated  in  Chapter  3  are 
compiled  together.  This  work  evaluated  computational  requirement  for  various  DOA 
algorithms  using  common  data  and  assumptions.  Computational  requirements  for  these 
algorithms  are  presented  in  Table  4.1.  A  review  and  challenges  of  these  algorithms  are 
summarized  in  Table  4.2. 

These  wide-band  DOA  algorithms  use  following  computational  steps: 

•  Generation  of  wide-band  signals 

•  Conversion  of  time  domain  signals  into  frequency  domain  via  FFT. 

•  Computation  of  covariance  matrices  in  frequency  domain. 
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•  Computation  of  eigenvalues  and  eigenvectors 

•  Computation  of  initial  estimates  of  number  of  sources 

•  Computation  of  initial  estimates  of  DOA 

•  Computation  of  transformation  matrices  and  focusing  on  central  frequency 

•  Computation  of  eigenvalues  and  eigenvectors 

•  Computation  of  number  of  sources 

•  Computation  of  final  estimates  of  DOA 

In  this  work  we  have  identified  common  computational  steps  and  they  can  be 
implemented  in  hardware.  Most  of  these  algorithms  follow  similar  computational  steps 
with  some  variations.  These  variations  can  be  adopted  if  we  use  re-configurable  approach 
and  implement  these  algorithms  in  FPGAs.  Chapter  5  presents  hardware  implementation 
of  these  algorithms. 
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NAME  OF 

THE 

TECHNIQUE 

COHERENT 

SIGNAL 

SUBSPACE 

ROTATIONAL 

SIGNAL 

SUBSPACE 

BINLINE  AR 
TRANSFORMATION 

BEAMFORMING 

INVARIANCE 

STEERED 

COVARIANCE 

SPATIAL 

RESAMPLING 

TOPS-  DOA 
ESTIMATOR 

FDFIB- 

BEAMFORMER 

Authors 

Wang  & 

Kaveh  [6] 

Hung  &  Kaveh 
[7] 

Shaw  [8] 

Ta-Sung  Lee  [19] 

Krolik  & 

Swingler  [17] 

Krolik  & 
Swingler  [18] 

Yoon 

[20,38,40] 

Do-Hong,  Russer 
[21,23-24] 

Computational 

steps 

Compute  64 
point  FFT 

Compute  64 
point  FFT 

Compute  64  point  FFT 

Beamspace 
transformation 
into  J  frequency 
bins 

Separate 
frequencies 
using  FFT 

Insert  Kn  - 1 

zeros 

Perform 

interpolation 

FFT  for  each 
block  (J) 

Compute  FFT 

Compute  33 
Covariance 
matrices 

Compute  33 
Covariance 
matrices 

Compute  33 

Covariance  matrices 

Compute  weight 
using  Least 

Square  fit  method 

Cross  spectral 
density  matrix  K 
for  each 
frequency 

Convolve  with 
the  low  pass 

FIR  filter 

Compute  sensor 
output  for  pre¬ 
selected 
frequencies 

J  beamforming 

network 

operation 

initial  DOA 
estimate  using 
MUSIC 

initial  DOA 
estimate  using 
MUSIC 

J  Beamspace  data 
matrix  focused  on 
single  ref. 
frequency 

Steered  cov. 
matrix  R  for 
each  angle 

Perform 

decimation 

operation 

K  Covariance 
matrices 

Compute 

covariance 

matrix 

Computation 
of  Focusing 
matrices 

Computation  of 

Focusing 

matrices 

Computation  of 

Focusing  matrices 

Beamspace  data 
matrix  is  designed 
in  least  square 
sense 

Inverse  of 
steered 
covariance 
matrix  R  for 
each  angle 

Compute 
covariance 
matrix  for  B 
samples 

Compute 
eigenvalues  & 
eigenvectors  for 
each  frequency 

Eigen- 

decomposition 

Eigen- 

decomposition 

Eigen- 

decomposition 

Eigen- 

decomposition 

Eigen- 

decomposition 

Spatial  power 
spectral  estimate 

Z  for  each  angle 

Operations  are 
performed  for 
each  frequency 

Define  signal  & 
noise  subspace 
for  each 
frequency 

Number  of 

sources 

Number  of 

sources 

Number  of 

sources 

Number  of  sources 

Number  of 

sources 

Determine  peak 
positions  of  the 
power 

Eigen- 

decomposition 

Compute 
rotational  signal 
subspace 
focusing  matrix 

Perform  MUSIC 
on  selected 
frequency 

Compute 

DOA  using 
MUSIC 

Compute  DOA 
using  MUSIC 

Compute  DOA  using 
MUSIC 

Proposes  a  root 
MUSIC  algorithm 
for  DOA 

Number  of 

sources 

Eigen- 

decomposition 

Perform 

MUSIC  to 
compute  DOA 

Number  of 

sources 

Perform 

MUSIC  to 
compute  DOA 

Table  4.1:  Computational  requirements  of  various  wideband  DO  A  algorithms 
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NAME  OF 

THE 

TECHNIQUE 

TIME 

DOMAIN 

FIB  USING 

FIR 

CYCLOSTATIONARY 
BASED  COHERENT 
METHOD 

Authors 

Ward  [26,29- 
30] 

GelU  &  Izzo  [39] 

Computational 

steps 

Design  FIR 
filters 
(primary  & 
secondary) 

Compute  64  point  FFT 

Form  J 

beamforming 

networks 

Compute  Covariance 
matrices 

Covariance 

matrix 

Focusing  matrices 

Eigen- 

decomposition 

Find  weight  function 

Number  of 

sources 

R  matrix 

Perform 

MUSIC 

Singular  value 
decomposition 

Eigen- 

decomposition 

Number  of  sources 

Perform  MUSIC 

Table  4.1:  Computational  requirements  of  various  wideband  DO  A  algorithms  (continued  from  the  previous  page) 
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Tab 


Name 

Authors 

Review 

Challenges 

Coherent  signal 
subspace  method 

Wang  & 
Kaveh 

Utilizes  linear  transformation  using  focusing  matrices 

Coherent  sources  can  be  resolved,  superior  detection  and 
accuracy 

Vulnerable  to  source  location  bias 

Initial  DOA  estimates  are  required 

Rotational  signal 
subspace  method 

Hung  and 
Kaveh 

Uses  rotational  transformation  matrices 

Improves  bias  and  performance 

Additional  computation  steps  (SVD) 

Initial  DOA  estimates  are  required 

Coherent  interpolation 

Bienvenu 

Interpolates  the  wavefield  to  emulate  a  spatial  sampling  which  is 
adopted  to  individual  frequency 

It  is  applicable  to  linear  arrays  and  error  increases  with 
the  increase  in  the  spatial  frequency.  It  will  limit  the 
range  of  DOA. 

Efficient  Wide-band 
Source  Localization 
Using  Beamforming 
Invariance  Technique 

Lee 

Forms  J  beamforming  matrices 

Computes  weight  using  least  square.  Focusses  on  single 
frequency.  Uses  root  MUSIC 

Utilizes  parallel  algorithm  to  offset  heavy  polynomial 
computations. 

DOA  estimator  is  suboptimum 

Iterative 

Too  many  computations  and  transformations  from  one 
computational  step  to  another. 

Difficult  to  implement. 

Signal  subspace  DOA 
estimator 

Yoon, 

Kaplan, 

McClellan 

Assumes  that  the  J  frequency  bands  of  the  sources  are  known 
Computes  FFT  for  J  blocks,  compute  K  covariance  matrices 
Compute  eigenvalues  &  eigenvectors. 

Define  signal  &  noise  subspaces,  compute  focusing  matrices 
Compute  DOA 

Initial  DOA  estimates  are  required 

There  is  some  missing  information  in  this  work 

Cyclo  stationary 
based  coherent 
method 

GelU  & 
Izzo 

Uses  focusing,  transformation  matrices,  weight  function,  and 

SVD 

Restricts  to  cyclo  stationary  signals 

Spatial  Resampling 

Krolik  & 
Swingler 

Perform  operation  for  each  frequency  bin  (33) 

Insert  zeros  to  interpolate 

Low  pass  filter  using  convolution 

Decimate  signals 

Form  covariance  matrix 

Compute  DOA  with  MUSIC 

Design  of  interpolator  filter. 

Stability  of  various  parameters  and  their  selection 
Processing  is  done  for  each  frequency  bin 

Steered  Covariance 
matrices 

Krolik  & 
Swingler 

Performs  FFT  on  input  data.  Computes  covariance  matrix  for 
each  frequency.  Computes  steered  covariance  matrix  for  each 
angle.  Inverse  all  steered  covariance  matrices.  Finds  peak  of  the 
power.  No  eigenvalues  or  eigenvectors 

33  different  covariance  matrices 

90  different  steered  covariance  matrices  (2  matrix 
multiplication  for  each  angles) 

90  inverse  matrices 

Compute  maximum  power 

Extensive  computations  for  each  steered  angle 

FDFIB  Beamformer 

Do-Hong 
&  Russer 

Compute  FFT 

Form  J  beamforming  networks 

Compute  covariance  matrices 

Perform  MUSIC 

Requires  array  geometry  and  guessing  of  common 
frequency.  Computational  steps  are  not  very  clear. 

le  4.2:  Reviews  and  challenges  of  various  wideband  DOA  algorithms 
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Chapter  5 


Hardware  Implementation 

Various  wide-band  DOA  algorithms  available  in  the  literature  have  been  investigated. 
Comparative  studies  were  performed  from  the  computational  requirement  point  of  view.  The 
goal  of  this  study  to  select  an  algorithm  or  class  of  algorithm  which  will  be  suitable  for  further 
study  and  will  be  a  candidate  for  hardware  implementation  for  real  time  application.  It  is  clear 
from  this  study  that  there  are  common  computational  modules  in  these  algorithms  and  most  of 
them  use  following  computational  steps: 

•  Generation  of  wide-band  signals 

•  Conversion  of  time  domain  signals  into  frequency  domain  via  FFT. 

•  Computation  of  covariance  matrices  in  frequency  domain. 

•  Computation  of  eigenvalues  and  eigenvectors 

•  Computation  of  initial  estimates  of  number  of  sources 

•  Computation  of  initial  estimates  of  DOA 

•  Computation  of  transformation  matrices  and  focusing  on  central  frequency 

•  Computation  of  eigenvalues  and  eigenvectors 

•  Computation  of  number  of  sources 

•  Computation  of  final  estimates  of  DOA 

It  can  be  seen  that  these  computation  requirement  would  require  special  purpose  hardware  in 
order  that  it  would  get  executed  in  real  time.  There  would  also  be  a  need  to  exploit  parallel 
processing  to  speed  up  the  computation  process.  There  are  number  of  ways  to  build  hardware  for 
these  applications  and  some  of  the  options  are  described  as: 

•  Commercially  available  Digital  Signal  Processor 

•  Field  Programmable  Gate  Arrays  (FPGAs) 

•  Application  Specific  Integrated  Circuits  (ASIC) 
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Use  of  commercially  available  DSP  chips  is  a  viable  way  to  design  DSP  hardware.  These 
DSP  devices  can  either  be  programmed  in  C  or  Assembly  language.  DSPs  have  general  purpose 
instruction  set  which  are  mapped  to  the  architecture  in  an  optimal.  They  may  not  be  suitable  for  a 
specific  application.  Modern  DSPs  offer  on-chip  multiply  accumulate  unit,  multiple  memories, 
specialized  instruction  set  for  signal  processing  applications.  Generally  support  software  is 
available  from  the  DSP  manufacturer.  Some  times  there  is  access  to  design  libraries  and  also 
design  boards  are  also  available.  One  drawback  with  this  approach  is  that  the  algorithms  are 
executed  in  sequential  fashion  and  programs  are  stored  in  the  memory  along  with  the  data.  This 
could  create  a  bottleneck  in  terms  of  speed  of  execution.  DSPs  are  also  clocked  with  certain 
clock  frequency.  This  approach  provides  a  viable  way  to  have  a  proof  of  concept  hardware.  In 
order  to  gain  more  speed  parallel  processing  can  be  exploited  and  multiple  chips  can  be  used. 

Second  option  of  developing  application  specific  hardware  is  to  use  FPGAs.  They  contain 
thousands  of  look  up  tables  to  store  logic,  hundreds  of  I/O  blocks,  on-chip  memory,  on-chip 
multiplier  and  very  flexible  programmable  multi- standard  I/O  pins.  There  are  number  of 
manufacturer  namely  Actel,  Alterra,  Xilinx  and  others  who  provides  these  devices  along  with 
synthesis  and  design  tools.  Different  manufacturers  use  different  technologies  and  operational 
philosophies  for  their  devices.  They  provide  their  own  proprietary  design  tools.  The  choice  of 
use  of  a  particular  FPGA  device  family  and  manufacturer  depends  on  availability  of 
experimental  devices,  design  tools,  training  facilities,  personal  preferences,  application 
dependency,  customer  specification  and  most  of  all  available  expertise  in  an  organization.  There 
is  a  steep  learning  curve  for  these  devices.  Switching  design  from  one  family  to  another  family 
of  FPGA  devices  may  again  require  additional  training  and  gaining  experience. 
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FPGAs  are  capable  of  implementing  high  performance  DSP  algorithm  as  they  can  provide 
multiple  Giga  Operation  Per  Second  (GOPS).  FPGAs  are  flexible,  programmable  and  can  be  re¬ 
configured  infinite  number  of  times.  Another  advantages  of  using  FPGAs  is  that  they  can  be 
used  to  implement  true  parallel  processing.  A  parallel  algorithm  can  be  implemented  easily  in 
FPGAs.  This  can  help  to  gain  speed  and  cut  down  the  execution.  This  kind  of  facility  or 
flexibility  is  unavailable  in  commercially  available  DSPs  as  we  have  to  worry  about  scheduling 
of  the  parallel  tasks.  The  operating  system  should  also  be  capable  of  handling  parallel  processing 
operations  and  be  capable  of  handling  scheduling  and  load  balancing  to  truly  achieve  benefits  of 
parallel  processing.  This  is  not  a  problem  in  FPGAs  as  custom  parallel  algorithms  can  be 
configured  and  mapped  on  the  FPGAs.  Various  designs  containing  parallel  algorithms  can  be 
experimented  with  the  FPGAs  to  achieve  high  speed  operations  and  optimize  need  for  FPGA 
resources.  This  kind  of  approach  can  totally  avoid  need  of  schedulers  and  operating  system  with 
parallel  processing  capabilities.  There  FPGAs  offer  a  good  rapid  prototyping  option  for  computer 
intensive  DOA  estimation  applications. 

Third  option  of  developing  special  purpose  hardware  is  with  the  use  of  Application  Specific 
Integrated  Circuits  (ASIC).  This  approach  is  capable  of  producing  design  which  are  custom 
designed  for  speed  and  area.  They  could  be  hand  crafted  to  achieve  desired  goal.  Drawbacks  for 
this  approach  are  high  design  cost,  very  high  level  of  design  expertise,  and  one  time 
programmability.  Other  drawback  would  be  that  each  ASIC  would  then  need  to  be  handcrafted 
for  each  application.  Moreover  device  needs  to  be  manufactured  by  the  manufacturer  themselves. 
This  will  take  away  flexibility  of  user  programmability  and  ability  to  change  designs  at  the  last 
minute.  Moreover  this  option  is  a  very  expensive  one  compared  to  the  previous  one.  This 
approach  is  not  pursued  in  this  work. 
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Hardware  Design  using  Digital  Signal  Processors 

A  general  purpose  DSP  was  selected  as  an  appropriate  platform  for  implementation 
because  of  its  ease  of  programming.  Also,  the  DSP  is  the  best  suited  for  matrix  and  floating 
points  computations.  The  DSP  used  in  our  implementation  is  a  DIOPSIS™  740  by  Atmel. 
DIOPSIS™  740  (D740)  is  a  high  performance  dual-core  processing  platform  for  real  time 
applications.  The  D740  is  optimally  suited  for  floating  point  applications  complex  domain 
computations.  The  ARM7TDMI  embedded  microcontroller  core  is  equipped  with  several 
peripherals  and  on-chip  memories.  The  main  components  of  the  DSP  subsystem  are  the  core 
processor,  the  on-chip  memories  and  the  interfaces  to  and  from  the  ARM  subsystem.  The  mAgic 
DSP  has  four  on-chip  memory  blocks:  the  program  memory,  the  data  memory,  the  data  buffer, 
and  the  dual  ported  memory  shared  with  the  ARM  processor.  An  external  memory  interface 
multiplexes  the  data  accesses  and  the  program  accesses  to  and  from  the  external  memory.  The 
program  memory  stores  the  Very  Long  Instruction  Word  (VLIW)  program  to  be  executed  by 
mAgic. 

Multicore  Application  Development  Environment  (MADE)  is  an  Integrated  Development 
Environment  (IDE)  that  can  be  used  to  develop  D740  applications  [43-45].  It  includes  the  C 
compilers  for  both  ARM  and  mAgic  DSP  based  on  GNU  compiling  tools  named  as  GCC.  The 
magic  C  compiler  contains  a  DSP  library  composed  of  over  220  functions  such  as  Fast  Fourier 
Transform  and  HR  and  FIR  filter  creation.  The  JTST  board  [43-45]  is  low-cost,  stand-alone, 
general-purpose  module  that  provides  the  appropriate  resources  in  order  to  evaluate  D740  DSP 
performances  in  a  wide  range  of  applications.  The  JTST  board  provides  several  memories  and 
other  peripherals. 
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DIOPSIS™  740  (DSP)  was  used  as  an  implementation  vehicle  for  the  Coherent  Signal- 
Subspace  algorithm  proposed  by  Wang  &  Kaveh  [6].  This  technique  separates  the  wide 
frequency  band  into  narrowband  components.  The  data  set  for  this  algorithm  is  divided  into  64 
segments  and  each  segment  contains  64  samples.  It  also  assumes  a  uniform  linear  array  of  16 
sensors.  The  Coherent  Signal  Subspace  approach  proposed  by  Wang  &  Kaveh  [6]  follows 
following  computational  steps. 

•  Compute  64  sets  of  64-point  FFT 

•  Compute  33  Covariance  matrices  (16  by  16) 

•  Computation  of  initial  DOA  estimate  using  MUSIC  algorithm  [14] 

•  Computation  of  Focus  matrix 

•  Computation  of  number  of  sources 

•  Separation  of  Signal  &  Noise  subspaces 

•  Compute  DOA  using  MUSIC  algorithm 

The  covariance  matrices  for  all  the  frequency  components  are  estimated  and  combined  to 
form  a  single  focused  covariance  matrix.  The  narrowband  MUSIC  algorithm  can  then  be  applied 
to  the  resulting  focused  covariance  matrix.  It  can  be  seen  that  the  CSS  algorithm  is  based  on 
matrix  computations  and  orthogonal  transformations,  which  are  computationally  intensive.  The 
eigendecomposition  problem  is  a  very  important  part  of  this  DOA  estimation  algorithm.  Finding 
the  eigenvalues  and  eigenvectors  of  the  covariance  matrix  is  needed  to  construct  the  signal  and 
noise  subspaces  that  the  CSS  algorithm  will  use.  The  Householder  and  QR  algorithms  [43]  can 
be  used  to  compute  the  eigenvalues  and  eigenvectors  of  the  symmetric  covariance  matrix.  The 
Householder  algorithm  is  used  to  reduce  the  bandwidth  of  the  covariance  matrix  by  transforming 
it  into  tridiagonal  form.  The  eigenvalues  and  eigenvectors  can  then  be  computed  using  the  QR 
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algorithm.  The  computational  cost  needed  for  computing  the  eigenvalues  and  eigenvectors  of  a 
tridiagonal  matrix  will  be  much  smaller  than  that  of  the  original  symmetric  matrix. 

A  parallel  architecture  capable  of  improving  the  performance  of  the  coherent  signal 
subspace  algorithm  was  proposed.  The  computation  time  of  the  proposed  architecture  was 
measured  and  compared  with  the  single  DSP  implementation.  The  results  showed  that  the 
parallel  architecture  yielded  the  same  results,  while  providing  superior  performance.  One  of  the 
limitations  of  the  proposed  architecture  was  the  use  of  static  matrices.  This  implied  that  the 
number  of  sensors  and  sources  was  known  in  advance.  It  also  implied  that  the  system  would  not 
be  able  to  detect  a  greater  number  of  sources  without  major  modifications  in  the  source  code. 
Using  dynamic  matrices  would  allow  the  system  to  easily  adapt  to  the  number  of  sources  to  be 


detected.  Table  5.1  shows  the  comparison  between  the  computation  time  of  the  single  DSP  and 


the  computation  time  of  the  parallel  architecture. 


Task 

Single  DSP 

Parallel  architecture 

Cycles 

Seconds 

Cycles 

Seconds 

Covariance  matrix 

1400000 

0.014 

42000 

0.00042 

Householder 

2600000 

0.026 

30000 

0.0003 

QR 

100000000 

1 

2100000 

0.021 

Power  spectrum 

1400000 

0.014 

47000 

0.00047 

Total 

210000000 

2.1 

5300000 

0.053 

Table  5.1  :  Performance  results  for  single  DSP  and  parallel  architecture 


Details  of  DSP  based  implementation  are  provided  in  the  second  volume  of  this  report.  A 
conference  paper  based  on  this  design  has  been  accepted  for  presentation  at  the  2007  IEEE  Radar 


Conference.  A  copy  of  the  paper  is  provided  in  Appendix  I. 


FPGA  based  Design 
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FPGAs  can  be  used  to  implement  wide-band  DOA  algorithms.  This  approach  will  provide 
shorter  time  to  design,  chance  to  exploit  parallel  processing  to  squeeze  all  the  processing  power 
and  a  chance  to  reconfigure  to  accommodate  other  designs.  In  this  work,  we  have  chosen  to  use 
Xilinx’s  FPGAs  for  multiple  reasons.  Most  of  all,  it  was  easy  to  access  design  tools  to  work  with. 
A  Field  Programmable  Gate  Array  (FPGA)  is  a  general  purpose  integrated  circuit.  It  is  user 
programmable  and  provides  flexibility  and  reconfigurability.  Xilinx’s  FPGAs  [46]  use  static 
RAM  to  keep  their  configuration  environment.  A  bitstream  file  is  created  and  downloaded  into 
RAM  for  FPGAs.  They  are  used  to  provide  their  configuration.  If  another  design  need  to  be  used 
then  another  file  can  be  downloaded  as  bitstream  file  and  FPGA  will  use  that  configuration. 

Xilinix  provide  excellent  design  tools  and  these  tools  are  also  available  for  research  work. 
Main  design  tool  consists  of  behavioral  simulation,  synthesizer  to  perform  functional  simulation, 
implementation  for  timing  simulation  and  download  feature  for  in-circuit  verification.  There  are 
two  mechanisms  to  provide  initial  design  specifications  to  the  design  tool. 

First  of  all  we  need  to  describe  any  design  in  a  hardware  descriptive  language.  In  our  case  we 
use  Very  High  Speed  Integrated  Circuit  Hardware  Descriptive  Language  (VHDL)  which  is  an 
industry  and  IEEE  standard.  The  VHDL  code  is  fed  to  Xilinx’s  synthesizer  which  performs 
synthesis  and  provides  mapping,  timing,  placement  in  a  bit  stream  file.  Xilinx  tool  [46]  will 
transform  VHDL  codes  into  standard  FPGA  components  such  as  Look  Up  Tables  (LUTs), 
outputs  F  and  Gs  from  the  LUTs  and  FFs.  It  would  yield  a  fuse  file  to  program  the  FPGA. 
Synthesis  tools  should  be  also  be  fast,  cost  effective  and  technology  independent.  This  approach 
reduces  risk  and  last  minute  changes  can  be  made  to  the  design.  It  also  optimizes  for  area  or 
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speed.  If  one  is  satisfied  with  the  synthesis  then  one  can  move  on  to  downloading  and  perform 
hardware  development  using  actual  FPGA  chip. 

One  drawback  of  this  approach  is  that  the  Xilinx  synthesizer  does  not  recognize  floating 
point  and  complex  numbers.  Wide-band  DOA  algorithms  heavily  rely  on  floating  point  and 
complex  numbers.  It  would  require  translating  of  all  these  data  types  to  integer  using  other 
software  available  else  where.  The  problem  with  this  approach  is  that  we  can  quickly  loose 
control  of  our  design  and  get  bog  down  with  so  many  details.  The  other  problem  is  that  this 
approach  requires  very  high  skill  level  of  VHDL  programming.  It  is  not  advisable  to  pursue  this 
approach  due  to  inefficiencies  in  the  synthesizer. 

There  is  an  alternate  mechanism  which  is  the  result  of  partnership  between  Mathworks  [47] 
and  Xilinx  [46]  Corporations.  They  provide  a  System  Generator  software.  This  software  would 
take  a  design  in  MATLAB/Simulink  code  and  system  generator  would  provide  a  VHDL  code 
which  that  can  be  used  with  the  Xilinx’ s  design  tools  (synthesizer).  This  will  by  pass  writing  of 
VHDL  code.  System  generator  also  provides  converter  blocks  for  converting  floating  point 
numbers  to  standard  bit-vectors  and  vice  versa.  This  feature  eliminates  the  need  for  special 
routines  for  dealing  and  redefining  the  floating  point  and  complex  numbers. 


System  Generator  [47]  is  used  to  create  design  in  Simulink,  simulation  capabilities  in  a 
graphical  environment.  It  allows  connection  of  subsystems  to  form  a  larger  system.  System 
Generator  has  a  library  of  blocks  some  times  called  cores  and  allows  incorporation  of  user 
defined  blocks  in  VHDL.  Library  provides  blocks  for  multipliers,  arithmetic  operations. 
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memories,  buffers,  up/down  data  rate  converters,  counters,  input/output  ports  and  other  blocks 
etc.  It  provides  a  capability  to  implement  the  designing  onto  FPGAs.  The  output  of  the  System 
Generator  is  VHDL  which  is  used  by  the  synthesizer. 

System  generator  approach  is  very  attractive  and  is  an  excellent  technique  for  designing 
specialized  hardware  for  DSP  application.  It  would  be  extremely  beneficial  for  us  to  convert  our 
wide-band  DOA  algorithm  code  and  re-develop  Simulink  code.  System  generator  also  offers 
some  pre-developed  cores  which  can  be  used  for  initial  rapid  prototyping.  This  work  is  using 
System  Generator  software  for  developing  FPGA  based  hardware.  The  hardware  design  work 
using  FPGAs  is  in  progress. 


FFT  Implementation 

First  computational  block  in  most  of  the  wide-band  DOA  algorithms  is  computation  of  Fast 
Fourier  Transform  (FFT).  Most  of  algorithms  and  our  MATLAB  program  use  64-point  FFT. 

This  FFT  block  is  used  continuously.  It  would  be  appropriate  to  investigate  best  way  to  design 
FFT  block.  Its  implementation  is  well  researched,  documented  and  widely  available  in  the 
literature.  There  are  number  of  DSP  and  other  chips  are  also  available  for  its  computation.  Most 
of  these  high-performance  FFT  chips  employ  parallel  arithmetic  units  and  cascaded  structures 
with  varying  processing  time.  Cascaded  structures  provide  pipelining  and  parallel  processing 
capabilities  and  provide  good  cost-performance  tradeoff.  First  of  all  it  is  proven  and  documented 
that  FFT  computation  using  radix-4  would  provide  efficient  implementation  when  compared  to 
radix-2  implementation.  It  is  proposed  that  we  should  use  radix-4  implementation  which  would 
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be  appropriate  in  our  case  as  the  only  restriction  in  radix-4  algorithm  that  the  number  of  points 


should  be  power  of  4.  It  is  true  in  our  case  as  we  would  like  to  compute  64-point  FFT  and  it  is 
power  of  4.  Literature  search  was  conducted  and  two  following  interesting  work  were  found: 

•  COBRA:  A  100-MOPS  Single-Chip  Programmable  and  Expandable  FFT  by  Tom  Chen, 
Glen  Sunada,  and  Jian  Jin  [27] 

•  An  Expandable  Column  FFT  Architecture  Using  Circuit  Switching  Networks  by  Tom 
Chen  and  Li  Zhu  [28] 

These  techniques  offer  interesting  structure  that  need  to  be  modified  and  adopted  for  FPGA 
implementation.  Following  building  blocks  would  be  required. 


1 .  An  array  of  16  bit  radix-4  butterfly  processors. 

2.  128*128  crossbar  switch  matrix 

3.  128-element  data  exchange  block 

4.  128-element  input/output  (I/O)  memory 

5.  Controller 


This  approach  would  require  design  of  highly  parallel  butterfly  processor  and  an 
interconnection  mechanism  using  crossbar  switch  matrix  concept.  The  I/O  memory  block  should 
be  divided  into  two  sub-blocks:  input  memory  block  and  the  output  memory  block.  The  memory 
sub-block  is  further  subdivided  into  two  sections:  one  for  real  part  of  a  64-point  complex  input 
vector  and  another  for  the  imaginary  part.  Loading  of  input  data,  circulating  of  intermediate  data 
and  sending  out  the  output  data  should  be  considered  and  design  in  an  optimal  way.  This  detailed 
design  work  has  not  been  pursued  and  at  this  time  we  will  use  available  block  from  the  System 
Generator  for  FFT  computation. 

The  work  on  other  computational  block  is  in  progress  and  would  be  reported  at  the  end  of 
this  contract. 


103 


Chapter  6 


Conclusions 

This  work  performed  a  review  of  wide-band  DOA  algorithms  in  the  literature  which  was 
accumulated  for  more  than  30  years  [1-41].  There  are  more  than  fifty  publications  for  wide¬ 
band  detection  of  Direction  of  Arrival  (DOA)  algorithms  which  are  available  in  the  literature. 
We  have  reviewed  the  most  relevant  one  and  have  not  reviewed  others  which  will  not  be 
applicable  or  suitable  from  hardware  implementation.  These  algorithms  were  generally  presented 
in  a  very  complex  or  condense  form,  which  are  not  easily  understandable  for  people  who  are 
outside  that  narrow  field.  One  of  the  reason  for  their  complex  representation  is  due  to  their 
publication  in  IEEE  transactions  and  conference.  These  transactions  generally  prefer  highly 
mathematical  papers  and  sometimes  authors  insert  mathematics  so  chances  of  their  papers  are 
increased.  Another  reason  for  condense  reporting  is  that  these  papers  face  a  page  limit.  Therefore 
algorithms  need  to  be  accommodated  within  those  guidelines  and  also  comply  with  the  reviewers 
comments.  One  unfortunate  thing  happens  in  this  process  that  essential  information  does  not  get 
into  the  papers  and  there  is  always  missing  information.  This  missing  information  is  acute  in  our 
case  as  we  are  looking  for  hardware  implementation  point  of  view  and  we  are  ignoring  details  of 
statistical  results  and  errors  which  are  irrelevant  in  our  case.  We  are  willing  to  sacrifice  small 
amount  of  error  in  order  to  accomplish  the  goal  of  implementing  them  in  hardware  for  real  time 
applications. 
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We  have  discovered  a  class  of  computational  requirements  that  would  be  required  in  all 
these  algorithms  as  was  summarized  in  Table  4.1.  We  have  also  given  reviews  and  challenges  in 
Table  4.2. 

We  were  able  to  cut  through  all  the  mathematics  and  converted  algorithms  into  simple 
arithmetic  operations.  This  step  is  very  useful  in  visualizing  an  architecture.  We  have  filled  a  gap 
between  the  design  of  computer  hardware  especially  special  purpose  parallel  architectures  and 
available  algorithms  for  various  interdisciplinary  problems. 

This  work  was  the  first  step  in  sorting  out  which  algorithm  is  appropriate  for  further 
study  and  for  its  hardware  implementation  for  real  time  applications.  We  have  developed  some 
hardware  as  described  in  Chapter  5  and  Volume  2  of  this  report.  Work  is  in  progress  for 
implementation  of  identified  computational  steps. 

This  work  can  be  extended  to  develop  re-configurable  test-bed  environment  for 
investigative  studies  for  various  algorithm.  The  re-configurable  test-bed  would  be  useful  to  study 
timing,  memory,  hardware  requirement  and  accuracy  of  results  from  various  algorithms.  This 
test-bed  would  also  be  useful  in  evaluating  different  number  of  sensors  and  different  kind  of 
sensors.  This  test-bed  would  also  be  a  technology  scalable  system  and  would  become  useful  in 
deployment  hardware. 
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Abstract 

There  is  a  need  for  the  computation  of 
Direction  of  Arrival  (DOA)  for  wideband  sources  for 
number  of  applications.  There  are  number  of 
algorithms  available  in  the  literature  for  wideband 
case  however  we  focus  on  Coherent  Signal  Subspace 
method  proposed  by  Wang  and  Kaveh.  Most 
algorithms  available  in  the  literature  follow  some 
variation  of  CSS  algorithm  thereby  increasing 
computational  complexity  in  an  effort  to  get  more 
accurate,  statistically  stable  and  unbiased  DOA.  We 
are  interested  in  computing  DOA  for  wideband 
sources  in  real  time.  We  chose  CSS  algorithm  and 
investigated  possibility  of  implementing  it  using 
commercially  available  Digital  Signal  Processor 
(DSP)  in  an  effort  to  achieve  real  time  capability. 
DSPs  offer  flexibility,  ease  of  development  of 
embedded  system,  reduces  design  cost  and  offers  use 
of  high-level  programming  language  such  as  C.  In 
this  work,  we  propose  a  DSP  based  architecture  for 
detecting  and  estimating  the  DOA  of  wideband 
sources.  It  is  known  that  DOA  algorithms  require 
computation  of  eigenvalues  and  eigenvectors.  It 
would  be  best  to  find  computational  friendly 
algorithm  for  the  computation  of  eigenvalues  and 
eigenvectors.  In  this  work  eigenvalues  and 
eigenvectors  are  computed  using  well  known 
Householder  and  QR  algorithms.  CSS  algorithm  is 
then  implemented  in  C  and  executed  on  DIOPSIS™ 
740  by  Atmel.  DIOPSIS™  740  (D740)  is  a  high 
performance  dual-core  processing  platform  for  real 
time  applications.  The  CSS  algorithm  was  then 
parallelized  and  a  parallel  architecture  was  then 
developed.  This  paper  presents  parallel  architecture 
using  DIOPSIS™  740  (D740)  and  computes 
performance  parameters. 


I.  Introduction 

Array  processing  has  been  an  important  part 
of  signal  processing  in  the  past  few  decades  [1-2]. 
The  array  consists  of  sensors  located  at  different 
points  in  space  with  respect  to  a  reference  point. 
Direction  of  Arrival  (DOA)  denotes  the  direction 
from  which  the  wave  fields  arrive  at  the  sensor  array. 
The  goal  in  DOA  detection  and  estimation  is  to 
accurately  determine  the  number  of  sources 
producing  waveforms  and  the  locations  of  those 
sources.  There  are  number  of  publications  available 
in  the  literature  for  detection  of  directional  of  arrival 
for  wideband  sources.  [1-15].  Wang  and  Kaveh  [1] 
proposed  Coherent  Signal-Subspace  (CSS)  method 
for  detection  of  DOA  for  wideband  sources.  This 
technique  separates  the  wide  frequency  band  into 
narrowband  components.  The  data  set  for  this 
algorithm  is  divided  into  64  segments  and  each 
segment  contains  64  samples.  We  also  assume  a 
uniform  linear  array  of  16  sensors.  The  Coherent 
Signal  Subspace  approach  proposed  by  Wang  & 
Kaveh  [1]  follows  following  computational  steps. 

9.  Compute  64  sets  of  64-point  FFT 

10.  Compute  33  Covariance  matrices  (16  by  16) 

11.  Computation  of  initial  DOA  estimate  using 

MUSIC  algorithm  [14] 

12.  Computation  of  Focus  matrix 

13.  Computation  of  number  of  sources 

14.  Separation  of  Signal  &  Noise  subspaces 

15.  Compute  DOA  using  MUSIC  algorithm 

The  covariance  matrices  for  all  the 
frequency  components  are  estimated  and  combined  to 
form  a  single  focused  covariance  matrix.  The 
narrowband  MUSIC  algorithm  can  then  be  applied  to 
the  resulting  focused  covariance  matrix.  It  can  be 
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seen  that  the  CSS  algorithm  is  based  on  matrix 
computations  and  orthogonal  transformations,  which 
are  computationally  intensive.  The 
eigendecomposition  problem  is  a  very  important  part 
of  this  DOA  estimation  algorithm.  Finding  the 
eigenvalues  and  eigenvectors  of  the  covariance 
matrix  is  needed  to  construct  the  signal  and  noise 
subspaces  that  the  CSS  algorithm  will  use.  The 
Householder  and  QR  algorithms  [16-17]  can  be  used 
to  compute  the  eigenvalues  and  eigenvectors  of  the 
symmetric  covariance  matrix.  The  Householder 
algorithm  is  used  to  reduce  the  bandwidth  of  the 
covariance  matrix  by  transforming  it  into  tridiagonal 
form.  The  eigenvalues  and  eigenvectors  can  then  be 
computed  using  the  QR  algorithm.  The 
computational  cost  needed  for  computing  the 
eigenvalues  and  eigenvectors  of  a  tridiagonal  matrix 
will  be  much  smaller  than  that  of  the  original 
symmetric  matrix. 

II.  DSP  Implementation 

To  implement  this  algorithm  in  real-time, 
hardware  capable  of  executing  millions  of  operations 
per  second  is  required.  A  general  purpose  DSP  was 
selected  as  an  appropriate  platform  for 
implementation  because  of  its  ease  of  programming. 
Also,  the  DSP  is  suitable  for  matrix  and  floating 
point  computations.  The  DSP  used  in  our 
implementation  is  a  DIOPSIS™  740  by  Atmel. 
DIOPSIS™  740  (D740)  is  a  high  performance  dual¬ 
core  processing  platform  for  real  time  applications 
[18].  The  D740  is  optimally  suited  for  floating  point 
applications  complex  domain  computations.  The 
ARM7TDMI  embedded  microcontroller  core  is 
equipped  with  several  peripherals  and  on-chip 
memories.  The  main  components  of  the  DSP 
subsystem  are  the  core  processor,  the  on-chip 
memories  and  the  interfaces  to  and  from  the  ARM 
subsystem.  The  mAgic  DSP  has  four  on-chip 
memory  blocks:  the  program  memory,  the  data 
memory,  the  data  buffer,  and  the  dual  ported  memory 
shared  with  the  ARM  processor.  An  external  memory 
interface  multiplexes  the  data  accesses  and  the 
program  accesses  to  and  from  the  external  memory. 
The  program  memory  stores  the  Very  Long 
Instruction  Word  (VLIW)  program  to  be  executed  by 
mAgic. 

Multicore  Application  Development 
Environment  (MADE)  is  an  Integrated  Development 
Environment  (IDE)  that  can  be  used  to  develop  D740 
applications  [18].  It  includes  the  C  compilers  for  both 
ARM  and  mAgic  DSP  based  on  GNU  compiling 
tools  named  as  GCC.  The  magic  C  compiler  contains 
a  DSP  library  composed  of  over  220  functions  such 


as  Fast  Fourier  Transform,  HR  and  FIR  filter 
creation.  The  JTST  board  [18]  is  low-cost,  stand¬ 
alone,  general-purpose  module  that  provides  the 
appropriate  resources  in  order  to  evaluate  D740  DSP 
performances  in  a  wide  range  of  applications.  The 
JTST  board  provides  several  memories  and  other 
peripherals. 

CSS  computation 

There  are  33  covariance  matrices  in  this 
algorithm  and  the  covariance  matrix  for  each 
frequency  component  is  computed.  The  main  issue 
resides  in  forming  the  matrix  X  containing  all 
samples  of  one  frequency  components.  Samples 
related  to  each  frequency  are  spread  across  the  64 
matrices  and  need  to  be  put  in  the  same  matrix.  This 
can  be  done  by  using  the  memory  transfer  function  to 
transfer  each  of  the  data  matrices  and  then  store  the 
row  corresponding  to  the  frequency  needed.  Figure  1 
shows  that  computation  process  for  the  frequency 

component  W0  . 

Using  the  periodicity  and  symmetric 
properties  of  the  FFT,  it  is  possible  to  have  all  the 

information  needed  by  selecting  frequencies  W0 

to  w32  .  The  process  above  is  repeated  for  each  of  the 

33  covariance  matrices.  The  eigenvalues  and 
eigenvectors  of  a  symmetric  matrix  are  computed 
using  the  Householder  and  QR  algorithms  [16-17]. 
However,  these  algorithms  cannot  be  directly  applied 
to  the  covariance  matrix  as  it  is  a  Hermitian  matrix. 
Solving  this  problem  requires  the  conversion  of  the 
Hermitian  matrix  into  a  real  symmetric  matrix.  It 
should  be  noted  that  Hermitian  matrices  have  real 
eigenvalues.  The  power  spectrum  is  computed  to 
obtain  an  initial  estimate  for  the  DOA.  The  column 
vectors  forming  the  noise  matrix  are  the  eigenvectors 
associated  with  the  M-D  lowest  eigenvalues,  where 
M  is  the  number  of  sensors  and  D  is  the  number  of 
sources.  The  number  of  sources  has  been  determined 
using  the  MDL  algorithm  and  initial  DOA  estimates 
are  obtained  [14-15].  The  computational  process 
continues  for  the  computation  of  focus  matrix, 
computation  of  eigenvalues  &  eigenvectors, 
computation  of  number  of  sources  and  then  final 
computation  of  DOA  estimate. 

III.  Parallel  Architecture  for  Coherent  Signal  Subspace 
Algorithm 

A  method  for  computing  the  covariance 
matrices  and  other  modules  using  a  single  DSP  was 
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Figure  1 :  Data  transfer  scheme  from  external  memory  to  local  memories 


Figure  2:  Parallel  architecture  for  coherent  signal  subspace  algorithm 
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presented  earlier.  This  method  can  be  improved 
upon,  by  using  one  DSP  to  perform  computation 
related  to  each  frequency  bins.  This  would 
require  computation  of  33  covariance  matrices. 
All  33  DSPs  share  the  same  external  memory. 
Figure  2  shows  its  parallel  architecture. 
Processors  can  initiate  a  DMA  transfer  from  the 
external  memory  to  their  local  memory.  The 
address  used  during  the  transfer  is  the  address  of 
the  row  containing  the  samples  at  a  particular 
frequency  and  the  number  of  elements 
transferred  is  the  number  of  elements  contained 
in  the  row.  The  DMA  transfer  for  the  first 
matrix  is  shown  above  in  Figure  3. 

After  a  DMA  transfer  is  complete, 
another  transfer  will  be  initiated  for  a  different 
segment  matrix.  This  process  is  repeated  for  all 
64  matrices  (segments)  in  external  memory.  The 
execution  time  will  be  the  same  for  each 
processor.  This  execution  time  can  be  further 
reduced  by  computing  only  the  lower  triangular 
part  of  the  covariance  matrix.  Since  the 
covariance  matrix  is  Hermitian,  the  information 
contained  in  the  lower  triangular  part  is  sufficient 
to  form  the  entire  matrix.  Using  a  single  DSP 
approach,  a  total  of  64  x  33  DMA  transfers  of 
256  elements  each  were  required.  The 
parallelized  process  only  requires  64  transfers  of 
16  elements  each,  which  should  result  in 
increased  performance. 

Using  certain  properties  of  the 

Householder  matrices,  it  is  possible  to  parallelize 
the  process.  This  observation  will  allow  us  to 
create  parallel  architectures  for  both  the 

Householder  and  QR  methods.  The  Householder 
method  consists  of  transforming  a  symmetric 
matrix  into  tridiagonal  form  using  the  following 
transformations 

b  =  h„_2...h2h1ah1h2  h„_2 

The  orthogonal  transformation  will  be 

accumulated  in  a  matrix  Q  in  order  to  recover  the 
eigenvectors  of  the  original  matrix  A.  This  can 
then  be  written  as  B  =  QAQ .  It  can  be  shown 

that  the  product  HA  could  be  computed  using  n 
parallel  processors  if  H  was  a  Householder 

matrix.  As  a  result,  E^F^A  and 

H„_2  •  •  •  H2Hj  A  can  also  be  computed  with  n 

processors.  The  tridiagonal  matrix  will  be  the 
result  of  the  product  of  two  sequences.  After 
execution  of  the  Householder  algorithm,  matrix 


A  is  the  tridiagonal  matrix  and  matrix  Q  contains 
the  product  of  all  the  orthogonal  transformations. 

The  QR  method  is  based  of  the  use  of 
the  following  orthogonal  transformations 

R  =  Hn_j . . .  H2HjA  and 
Q  -  H„_j . . .  H2Hj 

We  can  compute  R  and  Q  using  n  parallel 
processors.  The  next  step  involves  the 
computation  of  the  matrix  A  =  RQ  and  the 
computation  of  a  matrix  X  containing  the 
product  of  the  orthogonal  transformations 

Q1Q2  *  **Qm’w^ere  m  number  of 

iterations.  These  computations  can  be  done  using 
a  single  DSP.  This  process  will  start  over  until 
the  number  of  iterations  required  for  a  good 
approximation  has  been  reached.  After  execution 
of  the  QR  algorithm,  the  matrix  A  contains  the 
eigenvalues  of  the  original  matrix  along  its 
diagonal  and  matrix  X  contains  the  eigenvectors 
of  the  tridiagonal  matrix.  The  eigenvectors  of 
the  original  matrix  can  be  obtained  be 
multiplying  the  matrix  Q  obtained  in  the 
Householder  process  and  matrix  X. 

The  power  spectrum  needs  to  be 
computed  for  every  angle  between  0  and  90 
degrees.  Each  spectrum  value  can  be  computed 
by  sending  the  matrix  containing  the 
eigenvectors,  the  number  of  sources  and  the 
angle  of  arrival  to  a  specific  processor.  The  DSP 
can  then  compute  the  power  spectrum  and  send 
the  results  back  to  external  memory.  By  sending 
3  angles  values  to  30  of  the  33  available  DSPs,  it 
is  possible  to  reduce  the  time  needed  to  compute 
the  power  spectrum  by  a  factor  of  30. 

VI.  Simulation 

In  order  to  demonstrate  the  performance 
of  the  DO  A  algorithm  for  wideband  signals,  a 
uniform  linear  array  of  sixteen  equally  spaced 
Omni-directional  sensors  was  used.  The  spacing 
c 

between  sensors  is  - ,  where  c  is  the  velocity 

Vo 

of  propagation  and  f 0  is  the  central  frequency. 

Two  wideband  sources  at  0X  and  02  were 
assumed.  The  signals  are  stationary  zero  mean 
band  pass  white  Gaussian  processes  with  central 
frequency  fQ  =  100 Hz  and  bandwidth  B  =  40 Hz  . 
The  array  noise  is  also  stationary  zero  mean  band 
pass  with  the  same  pass  band  as  the  signal  with  a 
SNR  of  lOdB  at  each  sensor.  Source  signals  and 
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the  noise  are  random  processes  with  a  bandwidth 
of  40  Hz.  The  sampling  frequency  is  chosen  to 
be  300  Hz.  The  signal  will  be  observed  over  a 

period  of  To  seconds  and  To  will  be  divided 

into  k  —  64  segments.  On  each  of  those 
segments,  the  array  output  along  with 
corresponding  noise  will  be  decomposed  into 
narrowband  components  using  a  fast  Fourier 
transform.  The  total  number  of  samples  taken  by 
each  sensor  will  be  4096.  Simulation  data  is 
similar  to  the  method  described  by  Wang  & 
Kaveh  [1]. 

The  time  domain  samples  will  be 
transformed  into  frequency  domain  by  applying 
a  64  point  FFT  to  each  of  the  64  segments.  Due 
to  the  memory  limitation  of  the  DSP,  the  output 
data  will  be  saved  in  the  external  memory  as 
sixty  four  64x16  matrices,  where  each  matrix 
represents  the  output  for  each  of  the  64 
segments.  This  is  done  using  a  memory  transfer 
from  local  to  external  memory. 

The  previous  simulation  created  for 
single  DSP  was  used  again  for  the  case  of  the 
parallel  DSP  architecture.  This  was  done  to 
measure  DO  A  and  compare  their  performances. 
The  parallel  architecture  was  simulated  using  a 
single  DSP  by  executing  the  parallel  processes 
sequentially.  However,  the  performance  would 
be  measured  using  the  longest  running  process  in 
that  sequence  in  terms  of  number  of  clock  cycles. 
The  purpose  of  the  simulations  was  to  detect  and 

estimate  the  DOA  of  two  sources  located  at  20 

and  60  using  40  iterations  for  the  QR 
algorithm.  Figure  4  shows  the  simulations 
results.  It  can  be  seen  that  both  techniques 
yielded  the  same  results.  The  performance 
analysis  showed  that  the  most  computationally 
intensive  tasks  were  the  QR  algorithm,  the 
Householder  algorithm  and  the  power  spectrum 
computation.  Table  1  shows  the  comparison 
between  the  computation  time  of  the  single  DSP 
and  the  computation  time  of  the  parallel 
architecture. 

The  total  number  of  cycles  for  both  the  single 
DSP  and  the  parallel  architecture  includes  tasks 
that  could  not  be  parallelized  and  are  not  listed  in 
the  table.  The  total  execution  time  for  the  parallel 
architecture  is  0.05  seconds,  compared  to  2.1 
seconds  for  the  single  DSP.  This  is  due  to  the 
fact  that  the  execution  time  for  QR  algorithm, 
which  accounts  for  approximately  95%  of  the 


whole  process,  was  significantly  reduced  using 
the  parallel  architecture. 

VII.  Conclusion 

The  coherent  signal  subspace  algorithm 
was  chosen  for  the  implementation  because  of  its 
high  resolution;  the  platform  selected  was  an 
Atmel  Diopsis740  DSP.  A  parallel  architecture 
capable  of  improving  the  computational  time  of 
the  coherent  signal  subspace  algorithm  was 
proposed.  The  computation  time  of  the  proposed 
architecture  was  measured  and  compared  with 
the  single  DSP  implementation.  The  results 
showed  that  the  parallel  architecture  yielded  the 
same  results,  while  cutting  computational  time. 
One  of  the  limitations  of  the  proposed 
architecture  was  the  use  of  static  matrices.  This 
implied  that  the  number  of  sensors  and  sources 
was  known  in  advance.  It  also  implied  that  the 
system  would  not  be  able  to  detect  a  greater 
number  of  sources  without  major  modifications 
in  the  source  code.  Using  dynamic  matrices 
would  allow  the  system  to  easily  adapt  to  the 
number  of  sources  to  be  detected. 
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Figure  4:  Power  spectrum  for  final  DO  A  estimate  of  angles  0X  —  20  and#2  =  60  using  single  and 
parallel  DSP  approach _ _ _ 


Task 

Single  DSP 

Parallel  architecture 

Cycles 

Seconds 

Cycles 

Seconds 

Covariance  matrix 

1400000 

0.014 

42000 

0.00042 

Householder 

2600000 

0.026 

30000 

0.0003 

QR 

100000000 

1 

2100000 

0.021 

Power  spectrum 

1400000 

0.014 

47000 

0.00047 

Total 

210000000 

2.1 

5300000 

0.053 

Table  1  :  Performance  results  for  single  DSP  and  parallel  architecture 
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Abstract 

There  are  an  extensive  number  of  wide-band 
DOA  algorithms  available  in  the  literature.  If  one  needs  to 
implement  embedded  hardware  for  a  real-time  application, 
there  is  a  daunting  task  of  identifying  an  appropriate 
algorithm.  We  can  accept  a  small  amount  of  error  in  the 
hope  of  getting  hardware  which  can  execute  these 
algorithms  and  provide  results  in  real-time  for  a  given 
application.  This  work  is  a  first  step  in  sorting  out  which 
algorithm  will  be  appropriate  for  hardware  implementation 
and  presents  some  quantitative  comparisons  of  these 
algorithms. 

Introduction 

Passive  detection  of  objects  has  been  studied 
for  more  than  30  years  [1-13]  and  references  within 
them.  There  are  more  than  sixty  publications  for 
wide-band  detection  of  Direction  of  Arrival  (DOA) 
algorithms  which  are  available  in  the  literature. 
These  algorithms  are  generally  presented  in  a  very 
complex  or  condensed  form,  which  are  not  easily 
understood  by  people  who  are  outside  that  narrow 
field.  It  is  not  known  which  class  of  techniques 
would  be  appropriate  for  implementing  in  hardware 
and  would  be  useful  for  real-time  applications.  One 
has  to  cut  through  all  the  mathematics  and  convert 
algorithms  into  simple  arithmetic  operations  before 
the  architecture  can  be  visualized.  There  is  a  need  to 
bridge  a  gap  between  the  design  of  computer 
hardware,  especially  special  purpose  parallel 
architectures  and  available  algorithms  for  various 
interdisciplinary  problems. 

This  work  is  the  first  step  in  sorting  out 
which  algorithm  is  appropriate  for  further  study  and 
its  hardware  implementation  in  real-time 
applications.  Wide-band  DOA  algorithms  appropriate 
to  our  needs  have  been  reviewed  and  have  been 
simulated  in  MATLAB.  These  algorithms  can  be 
divided  into  following  categories: 

•  In-coherent  Signal  Subspace 

•  Coherent/rotational  Signal  Subspace  (CSS) 

•  CSS  using  other  focusing  matrices 


•  DOA  detection  using  beamforming 

•  Combination  of  beamforming  and  focusing 
approaches 

•  Use  of  ARM  A  model  and  Bayesian 
approaches 

•  Use  of  maximum  likelihood  algorithms 

•  Bilinear  Transformation  Method 

Signal  subspace  approaches  are  very  popular 
for  computing  DOA  for  both  narrow-band  and  wide¬ 
band  sources  [1-6].  One  problem  with  them  is  that 
they  may  not  be  optimal  but  they  produce 
computationally  efficient  algorithms.  Signal  subspace 
based  approaches  are  further  subdivided  into 
coherent  based  and  incoherent  based  signal  subspace 
approaches.  Incoherent  signal  subspace  approaches 
decompose  signals  into  individual  narrow-band 
frequencies  and  then  combine  them  to  produce  the 
final  results.  They  are  computationally  expensive. 
Some  of  the  subspace  approaches  require  initial 
estimates.  If  the  initial  estimate  is  not  accurate  then 
the  final  DOA  will  also  not  be  accurate  and  there  will 
be  issues  with  bias  and  variances. 

In  order  to  implement  a  DOA  algorithm  in 
hardware  for  real-time  applications,  it  is  important  to 
use  a  computationally  efficient  algorithm.  One 
approach  is  to  evaluate  the  computational 
requirements  of  currently  available  wide-band  DOA 
algorithms  and  select  one  of  them  for  hardware 
implementation. 

Wide-band  DOA  Algorithms 

Development  of  wide-band  DOA  algorithm 
started  with  the  incoherent  signal  subspace  approach 
which  is  a  brute  force  extension  of  the  narrow-band 
case  and  is  computationally  expensive.  Wang  & 
Kaveh  [5]  proposed  a  Coherent  Signal  Subspace 
Method  (CSSM)  which  creates  a  focusing  matrix 
using  a  single  frequency.  This  has  a  simple 
transformation  scheme  and  requires  an  initial 
estimate  of  the  DOA.  Results  are  reasonable  and 
produced  peaks  at  the  appropriate  DOA.  Hung  & 
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Kaveh  [6]  extended  their  previous  work  with  a 
promise  of  statistically  better  results.  They  introduce 
the  concept  of  a  Rotational  Signal  Subspace  (RSS) 
producing  a  more  complex  or  accurate  focusing 
matrix.  This  required  the  additional  computational 
step  of  singular  value  decomposition.  There  may  be 
little  improvement  in  accuracy  of  the  results.  Shaw 
[7]  extended  the  work  of  Wang/Hung  &  Kaveh  [5-6] 
and  introduced  a  bilinear  transformation  approach 
with  certain  approximation.  The  advantage  of  their 
work  is  that  it  does  not  require  an  initial  DOA 
estimate.  Their  algorithm  is  very  sensitive  to  certain 
assumptions  and  parameter  values  which  make  it 
unattractive  for  hardware  implementation. 

Krolik  et  al.  [8]  introduced  the  computation  of  a 
steered  covariance  matrix  for  each  angle,  then 
inverted  the  covariance  matrices  to  find  peaks  of  the 
power.  He  eliminated  the  need  of  eigen- 
decomposition  and  selection  of  an  initial  focusing 
angle.  This  technique  requires  computation  for  each 
steering  direction  6  increasing  the  computation 
requirements.  His  extended  work  [9]  uses 
interpolation  and  then  decimation  of  the  input  data 
followed  by  computation  of  the  covariance  matrix. 
The  technique  looks  good  but  some  work  needs  to  be 
done  to  resolve  some  of  the  parameters  as  their 
selection  may  be  tricky.  This  approach  may  not  be 
useful  in  a  generalized  case.  It  also  eliminates 
preliminary  DOA  estimates. 

Ta-Sung  Lee  [10]  decomposes  received  data 
using  bandpass  filters  into  J  frequency  beams.  His 
algorithm  performs  beamspace  transformation  and 
computes  weights  using  a  least  square  method.  It  then 
computes  a  beamspace  data  matrix  and  focuses  on  a 
single  reference  frequency  which  would  be 
something  similar  to  CSSM  method.  It  performs 
transformation  into  K  beamspaces  and  forms  the 
beamspace  data  matrix.  This  beamspace  data  matrix 
is  then  focused  on  a  single  reference  frequency  out  of 
J  frequency  bands.  The  design  of  the  beamspace  data 
matrix  is  described  which  first  requires  design  of 
beamforming  matrices.  The  weights  are  again 
designed  in  a  least  square  sense.  The  problem  is  then 
reduced  to  beamspace  data  correlation  and  noise 
matrices.  The  authors  then  apply  their  own  derived 
root  MUSIC  algorithm  which  could  be  substituted 
with  the  MUSIC  algorithm.  The  algorithm  does  not 
require  any  preliminary  DOA  estimates.  This  DOA 
estimator  is  suboptimum  according  to  the  author.  It 
would  also  result  in  degradation  in  estimation 
accuracy  at  low  SNR.  Even  with  the  true  DOA,  the 
DOA  estimates  may  not  be  exactly  identical  as  the 
signal  roots  may  not  lay  on  the  unit  circle.  This 
approach  provides  an  excellent  alternative  to  CSSM 


except  for  the  frequency  decomposition  and 
calculations  of  weight  and  beamspace  data  matrices. 


Tuan  Do-Hong  et  al.  [11]  first  compute  the 
FFT  then  form  beamforming  networks,  compute 
covariance  matrices,  and  then  perform  MUSIC.  They 
require  array  geometry  and  guessing  of  a  common 
frequency.  Daren  et  al.  [12]  performs  filter  and  sum 
beamforming  in  a  frequency  invariant  fashion.  They 
proposed  a  FIR  filter  design  and  then  computed 
covariance  matrices.  Three  FIB’s  were  designed 
using  FIR  filters  in  a  frequency  invariant  fashion 
covering  the  normalized  frequency  band  and  the 
spatial  sector.  This  approach  looks  very  attractive  and 
needs  some  more  work.  Filter  specifications  were  not 
provided.  We  need  to  find  a  way  to  design  these  FIR 
filter  coefficients  and  a  more  focused  way  to  compute 
them.  Sellone  [13]  proposed  a  new  way  of  designing 
focusing  matrices  which  has  the  same  robustness.  It 
does  not  need  initial  estimates  of  DOA.  He  also 
claims  that  the  computational  requirement  is  also 
reduced  when  compared  to  the  RSS  and  SST 
approach.  The  work  of  Yoon  [14]  requires  a 
preliminary  estimate  of  the  DOA  that  could  be  one  of 
the  angles  in  the  estimation  and  the  number  of 
sources  need  also  needs  to  be  estimated.  They 
compute  K  covariance  matrices,  eigenvalues  and 
eigenvectors.  They  form  signal  &  noise  subspaces 
and  compute  focusing  matrices.  Finally  they  find  the 
DOA.  They  have  also  named  this  technique  as  Test 
of  Orthogonality  of  Projected  Spaces  (TOPS). 


Algorithm 

Challenges 

CSSM  [5] 

Initial  DOA  estimates  are  required 

RSS  [6] 

Additional  SVD  step  and  initial  DOA  is 
required 

Beamforming 
Invariance  [10] 

DOA  estimator  is  suboptimum  and  iterative 

DOA  estimator 
[14] 

Initial  DOA  estimates  are  required 

Spatial 

Resampling  [9] 

Design  of  interpolator  filter. 

Stability  of  parameters  and  their  selection 
Processing  is  done  for  each  frequency  bin. 

Steered  Covariance 
matrices  [8] 

Steered  covariance  matrices  for  each  angle. 
Compute  maximum  power  for  each  steered 
angle 

FDFIB 

Beamformer  [11] 

Requires  array  geometry  and  guessing  of 
common  frequency. 

Table  2:  Challenges  of  wide-band  DOA  algorithms 


Computational  Requirements 

This  work  evaluated  the  computational 
requirement  for  various  DOA  algorithms  using 
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common  data  and  assumptions.  Computational 
requirements  for  these  algorithms  are  presented  in 
Table  1.  Their  implementation  was  also  explored. 
These  wide-band  DOA  algorithms  use  following 
computational  steps: 

•  Generation  of  wide-band  signals 

•  Conversion  of  time  domain  signals  into 
frequency  domain  via  FFT. 

•  Computation  of  covariance  matrices  in  the 
frequency  domain. 

•  Computation  of  eigenvalues  and 
eigenvectors 

•  Computation  of  initial  estimates  of  the 
number  of  sources 

•  Computation  of  initial  DOA  estimates 

•  Computation  of  transformation  matrices  and 
focusing  on  central  frequency 

•  Computation  of  eigenvalues  and 
eigenvectors 

•  Computation  of  number  of  sources  and  final 
estimates  of  DOA 

Conclusions 

In  this  work  we  have  identified  common 
computational  steps  and  they  can  be  implemented  in 
hardware.  Most  of  these  algorithms  follow  similar 
computational  steps  with  some  variations.  These 
variations  can  be  adopted  if  we  use  a  re-configurable 
approach  and  implement  these  algorithms  in  FPGAs. 
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Table  1:  Computational  requirements  of  various  wideband  DO  A  algorithms 
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